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FOREWORD 


This report describes the research performed during the period July 1990 - July 1991 under 
a joint research program between the Mechanical Engineering Department of the University 
of Maryland and the Buildings and Fire Research Laboratory of the National Institute of 
Standards and Technology. The research is conducted by Graduate Research Assistants of 
the ME Department under the joint supervision of Dr. di Marzo (UMCP) and Dr. Evans 


(CFR - NIST). This joint research program was initiated in January 1985. 


The long term objective of the study of droplet-solid interaction is to obtain information 
applicable to the extinguishment of fire through a droplet array (e.g. spray). The solids of 


concern include low thermal conductivity materials, typical of fire applications. 


Several important results were obtained in the first years of research. In particular, the 
modelling of the boundary condition at the liquid-vapor interface (at the droplet exposed 
surface) was validated with the data collected for water droplets evaporating on an 
aluminum block (diMarzo 1986a, 1986b, 1988). The simple model that describes the cooling 
effect induced in the aluminum by the evaporating droplet (diMarzo and Evans 1986, 
diMarzo et al. 1987) provided the input for the formulation of a preliminary multi-droplet 
model (diMarzo et aL 1987). An extensive experimental data base was also collected for 
water droplet evaporating on Macor (a glass like material able to withstand strong thermal 


stresses). An infrared thermographic technique has been developed to monitor the 


iv 


temperature distribution on the surface of Macor during the evaporative process. These 
data are used to validate the more complex model required to describe the coupled liquid- 


solid interaction for the Macor case. 


The research described hereafter constitutes portion of an extensive research program aimed 
at developing accurate droplet cooling models of burning solid fuel surfaces. Many studies 
have been performed to quantify the vaporization process of both single and multi droplet 
arrays impacting on hot surfaces. For the studies found in the published literature, the full 
span of the droplet vaporization processes are usually reported. These would include 
evaporation, nucleate boiling, film boiling and Leidenfrost transition. This present research 
is more limited in the span of vaporization process studied, being only concerned with the 
evaporation of single droplets on a hot surface. Limited experiments are also performed 


in the nucleate boiling range. 


This report consists of three different portions; each portion focuses on an independent issue 


of the research program. 


In Part A the major progresses in the numerical simulation of a single droplet evaporation 
are described. A comprehensive model for the prediction of evaporative cooling of solid 
surfaces induced by the impingement of single water droplets is presented. The model 
predicts the droplet evaporation and the solid surface cooling for materials with thermal 


conductivity spanning over more than two orders of magnitude. The model accurately 


predicts the total evaporation time and it is further validated with transient surface 
temperature measurements obtained by infrared thermography. The predictions are in 
excellent agreement with the experimental data. The spatial and temporal heat flux 
distribution under the evaporating droplet is studied. The extent of the droplet evaporative 
cooling is quantified by introducing a novel concept defining the droplet radius of influence. 
A closed fomnigoiatiad predicting the radius of influence is derived and tested against the 
experimental data and the model predictions. An empirical correlation for the prediction 
of the evaporation time is also presented. Insight into the evaporative cooling phenomena 
is provided for materials with various thermal conductivity. The relevant parameters are 


identified and their influence on the phenomena is assessed. 


In Part B the experimental program conducted with radiant heat input for single droplets 
evaporating on Macor is described. A new experimental set-up for the study of dropwise 
evaporation in a radiant heat transfer field has been designed, constructed and tested. The 
various issues of concern (steady state solid temperature distribution, radiant heater design 


and configuration, infrared background noise and post test data manipulation) are outlined. 


In Part C a novel model is presented, which extends the single droplet results to a multi- 
droplet system. Numerical results are obtained for different values of the impingement area, 
mass flow-rate of the spray and initial solid surface temperature. However, this beget of 
the research is still in progress and it will constitute an essential portion of the next grant 


period. 
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PART A 


SINGLE-DROPLET MODEL 


Al. INTRODUCTION 


The cooling of hot surfaces by droplet sprays has been studied extensively. Multi-droplet 
systems were investigated by Bonacina et al. (1979). He reinforced Toda’s (1972) conclusion 
that conduction in the liquid is the relevant heat transfer mode in absence of nucleate 
boiling. Convective heat transfer within the liquid droplet could be due to Rayleigh 
instabilities or more likely to surface tension. The Rayleigh instability driven motion should 
be ruled out due to the large response time which exceeds the droplet evaporation time. 
The surface tension induced convective motion is more plausible (see Ostrach and Pradhan, 
1978). However, the velocities that would results from this driving mechanism are very 
small (order of 10* meters-per-second). The temperature profiles in the liquid are fieene 
for such small internal motions. Flow visualization experiments were conducted by the 
authors with small ink tracers and no significant convective motion was detected even for 
large droplets under large thermal gradients. Sadhal and Plesset (1979) investigated the 
effect of solid properties. Grissom and Wierum (1981) further developed these ideas to 
define a range of conditions for spray evaporative cooling. Rizza (1981) formulated a 
thermal model for a semi-infinite solid subjected to droplet evaporation on its exposed 
surface. His main assumption of uniform and constant solid-liquid interfacial temperature 
is corroborated by Seki et al. (1978), and by the extensive experimental observations of 
Makino and Michiyoshi, see for example Makino and Michiyoshi (1979,1984). A simple 
model for single droplet evaporation was introduced by diMarzo and Evans (1989). This 


model is limited to evaporation at the liquid-vapor interface in complete suppression 
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of nucleate boiling. Photographic inspection of evaporating droplets (diMarzo and Evans, 
1987) revealed that the droplet shape can be approximated by a segment of a sphere. 
Zhang and Yang (1982) confirmed this assumption along with the more general behavior 
of droplet-surface interactions described by Simon and Hsu (1971) and by the extensive 


work of Wayner (see as an example Wayner, 1973). 


These studies focused mostly on water droplet deposited on high thermal conductivity 
materials. Thus, the assumption of uniform and constant solid-liquid interfacial temperature 
is reasonable and the modeling of the liquid region can be de-coupled from the treatment 
of the solid. This is not the case for dropwise evaporative cooling taking place on low 
conductivity materials. Abu-Zaid and Atreya (1989) measured the interfacial temperature 
at various locations and confirmed that the temperature changes during the transient. The 
thermal behavior of the surface during the evaporation, monitored by infrared thermography 
(Klassen et al., 1990), also reveals that a sharp change in the interfacial conditions is taking 


place in the later part of the droplet evaporation transient. 


The objective of this study is to provide a basis for the modeling of solid fuel fire 
extinguishment process and for the modeling of the cooling mechanism by sprinkler 
generated sprays in a fire environment. The thermal behavior of different materials 
subjected to dropwise evaporative cooling is studied. In particular, the effect of their 
different thermal conductivity (the range considered spans more than two orders of 


magnitude) is assessed. A model which describes the coupled thermal behavior of solid and 
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liquid during the evaporative transient has been formulated. This model is fully described 
by Tartarini et al. (1990) and will be the subject of a paper Sapa in preparation. The 
predictions of the model are presented and compared with the experimental results. Further 
analysis of numerical and experimental results is carried out to provide more insight into 
the evaporative cooling mechanism. A novel correlation for the prediction of the total: 
evaporation time is proposed. An approximate, closed-form solution is used to deduce the 


relevant normalized governing parameters. 


A2. PHENOMENOLOGY 


A droplet gently deposited on a hot surface wets a circular area of radius R which remains 
constant throughout all the evaporation. The water-solid contact angle decreases without 
reaching its minimum regressive value (see Simon and Hsu, 1971). Therefore, the shape is 
described by a segment of a sphere of constant base. The heat is conducted through the 
droplet from the solid surface to the liquid vapor interface where it evaporates the liquid 
and is convected and radiated into the ambient. The thermal conductivity of the solid 
determines the extent of the cooling and the magnitude of the local the temperature drop. 
As the evaporation progresses the thickness of the liquid layer, that constituted the droplet, 
decreases thus the heat transfer increases due to the reduced resistance. This fact implies 
that the heat flux at the droplet eee increases as time elapses and the portion of the liquid- 
vapor interface where most of the evaporation takes places extends gradually from the 
droplet edge to all the interfacial area as also the central portion of the droplet (at small 
radii) thins. Detailed description of phenomenology of the droplet evaporation and of the 
solid surface temperature behavior is provided by diMarzo and Evans (1987) and by Klassen 
et al. (1990). Experiments were conducted on aluminum and on Macor (a glass-like 
material). Some important physical properties of aluminum and Macor are reported in 
Table 1. The thermal conductivity of aluminum is about 180 W/m°C while Macor exhibits 
a thermal conductivity of 1.3 W/m°C. This difference of more than two orders of magnitude 
spans over most of the non-porous, solid materials range. The thermal conductivity more 


than the thermal diffusivity of the materials is the main parameter characterizing the 
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evaporative cooling behavior. 


Figure Al illustrates the differences between the aluminum and Macor behaviors. The 
temperature profiles on the surface of the two materials are compared (at different times 
during the transient t/r = 0.3 and t/r = 0.9) for an identical heat sink (30 ul water droplet 
evaporating in 95 seconds). Note that all the aluminum cooling is bounded within half 
degree Centigrade while the Macor cools over 30 degrees Centigrade. This huge difference 
in behavior explains why a constant and uniform interfacial temperature can be assumed for 
high conductivity materials and leads to meaningless results for a low thermal conductivity 


solid. 


The region of the solid surface affected by the droplet presence is confined to few droplet 
diameters for Macor while its effects extend on a large portion of the aluminum surface. 
The actual temperature drop at the solid surface is significant for Macor although it is very 
small for aluminum. These observations will be quantified in the following and provide 


insight on the cooling strategies to be used for different materials. 


The results obtained for high thermal conductivity solids are valid over a range of materials 
(copper, aluminum, steel). In all these cases, uniform and constant interfacial temperature 
can be assumed. De-coupling the liquid from the solid means that the liquid droplet 
behaves independently of the substrate as long as the interfacial temperature is close to the 


initial solid surface temperature. The interfacial temperature at the initial liquid-solid 
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contact can be obtained considering the exact solution available for two semi-infinite solids. 
Seki et al. (1978) bases his analysis on this consideration and suggests that the interfacial 
temperature should be: 


p= 21yP161%1 + TryPsCsks (A1) 
Cc 
¥P 161k; Ba ¥P3Csk, 


The relation between this interfacial temperature and the actual computed value will be 


discussed in the following and its use in the analysis will also be outlined. 


A3. THEORETICAL MODEL 


A3.1. Model description 
The modeling of the coupled solid and liquid thermal behavior is described by the transient 
conduction equation for both domains with the appropriate boundary conditions. The 


governing equations, with respect to the coordinate system depicted in Fig. A2, are: 


C] 


| 


ar = Beaver (A2) 
— =«,V*7 (A3) 


A linear temperature distribution in the solid is required as initial condition. The far-field 
temperature distribution must be unaffected by the droplet evaporation. Initially, the liquid 
droplet is in thermal equilibrium with the ambient. By introducing an overall heat transfer 
coefficient h at the exposed solid surface, the boundary conditions at the liquid-vapor 


interface and at the exposed solid surface can be written as: 


ala0rsel as Ry 270. 


(A4,A5,A6) 


Al ibin,. 2=.0: 


k, 7 = h(T,-T,) 


At the liquid-vapor interface a small portion of the heat conducted from below through the 
liquid is convected and radiated into the ambient. Most of the heat evaporates the liquid. 
To account for the evaporation, the vapor diffusion in the air is considered. The mass 
transfer coefficient is obtained with the Chilton-Colburn analogy (Chilton and Colburn, 
1934). Details of the derivation of this liquid-vapor interfacial condition are given by 
diMarzo and Evans (1987). Note that the overall heat transfer coefficient h account for all 
the heat transferred to the ambient while the mass transfer coefficient is related to the 
convective heat transfer h,,. The final formulation can be written as: 


-k, VT-A(T,-T,) = 
(A7) 


Extremely strong local thermal gradients at the drop initial contact anedindts the solution of 
this problem with conventional finite difference schemes. The solution is obtained by using 
a Boundary Element Method (BEM) for the solid region and a Control Volume Method 
(CVM) for the liquid region. The BEM is described in detail by Kavoosi et al. (1989). The 
adjoint equation to Eq. (A2) (see Wrobel and Brebbia, 1981) is satisfied by the following 


Green’s function: 


Aric ttct tne s eee ee 


(A8) 
_ (r-re)2+ (z- ze)? _ (r-re)2 + (z+z9)? 
e ance +e ance 


A linear combination of Eq. (A2) and of its adjoint equation, is integrated over the solid 
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and temporal domain. The transient terms are reduced to two spatial integrals at the initial 
time and at the present time. The integral at the initial time is eliminated with a proper 
choice of a transformed temperature as shown in Eq. (A9). The spatial integral at the 
present time is simply the integrand since the Green’s function reduces to a Dirac function. 


After defining: 


CZ OS ho Be (A9) 


3g 


and by using the Gauss theorem on the right hand side one obtains: 


(r-r*)? 
tre card chau 
BU ay te) S Ih f Vu(r*, t*) r* t#73/2 1, 222 | @ *%s* dregee 


4a,t* 
(A10) 

Details on the CVM are given by Tartarini et al. (1990). 

The governing equation for the liquid is cast in the following form: 


aT 4. ai 
LS dv=a, if fids (All) 


which is integrated for each elementary volume of the discretized liquid domain. These two 
equations, Eq. (A10) and Eq. (A11), are then combined in matricial form. By inverting the 
matrix, one obtains the transient temperature and heat flux distribution in the liquid and 
over the solid surface. The description of this complex numerical solution is the subject of 
a paper currently in preparation. 

The input to the code must prescribe: a) the droplet initial volume; b) the droplet shape 


parameter 8 as defined by Bonacina et al. (1979); c) the initial surface temperature; d) the 
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overall heat transfer coefficient and e) the convective heat transfer coefficient. The shape 
eae eeter and the two heat transfer coefficients are provided by correlations based on 
experimental data. The code predicts: a) the transient temperature profiles on the solid 
surface and at all liquid locations; b) the heat fluxes at all points of the liquid domain and 


on the solid surface; c) the total evaporation time and d) the transient liquid inventory. 


A3.2. Model validation 

The validation of the model predictions is obtained by comparing them with experimental 
results obtained by Kidssen et al. (1990). These experiments provide the temperature 
distribution during the evaporative transient of a water droplet deposited on Macor. Figure 
A3 illustrates Sine of these comparisons for various initial surface temperatures and droplet 
volumes at a time corresponding to 90 per cent of the total evaporation time. A comparison 
toward the end of the evaporation time is selected to maximize the discrepancies between 
predictions and data due to the parabolic nature of the problem. These local temperature 
profile comparisons (for r/R > 1) illustrate the typical model capabilities which are deemed 
satisfactory. A comprehensive comparison of the model with the data is obtained by looking 
at the overall predicted evaporation time for aluminum and Macor and at the experimental 
results obtained by Klassen et al. (1990) and by diMarzo and Evans (1989). Figure A4 
shows the calculated and measured evaporation time for both materials at various initial 
solid surface temperatures and droplet sizes. The agreement is remarkable in light of the 


fact that the experimental data are repeatable within ten per cent of the indicated value. 
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A4. RESULTS AND DISCUSSION 


A4.1. Approximate, closed-form solution 

In an effort to normalize the relevant parameters, an approximate, closed-form solution is 
sought. Carslaw and Jaeger (1959) show that, for a semi-infinite solid with a surface heat 
flux constant and uniform applied over a disk of radius R and no heat transfer on the 


remainder of the surface, the transient temperature profile on the solid surface is given by 
T= GE f, Jg(hr) J, (AR) erf(aja,e) (A12) 
Ss 


By introducing the normalized radius n = r/R and the normalized time 6 = (a, t)/*#/R and 

by setting q equal to the spatial and temporal average heat flux due to evaporation, one can . 

express the dimensionless temperature 8 as follows: 
Q = za £2) 8? B? = 


3 (A13,A14) 


4 fiona) J. (€) erf (8) = 


where f is the shape parameter defined by Bonacina et al. (1979) and JA is the Jakob 


number. These two parameters can be written as: 


B = r(S3)" (A15) 
ga = £hts7T) (A16) 
A 
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Figure A5 shows a comparison of the approximate solution with the data. The approximate 
solution under-predicts the temperature profile (for 7 > 1). However, the general trends 
are preserved and substantial advantage is gained with this normalization. This closed-form 


solution will also be used in the following to assess the droplet cooling effect. 


A4.2. Results 

The validated model is used to investigate the interfacial temperatures and heat fluxes under 
the droplet where direct measurements are unavailable. Figure A6 shows the liquid-solid 
interfacial fluxes for a typical case. It is important to note that the heat flux is not uniform 
nor constant during the evaporative process. The spatial distribution indicates that most of 
the evaporation takes place at the outer edge of the droplet. This results confirms the 
Endings of diMarzo and Evans (1989). Therefore it is mot surprising that the 
approximate solution, based on constant and uniform flux, fails to match the surface 
temperature profiles for n =1. 

Figure A7 illustrates the behavior of the spatial averaged heat flux for different materials. 
‘The onset of nucleate boiling temperature for these five cases is calculated by setting the 
contact temperature in Eq. (A1) to 100 °C. The resulting initial solid surface temperature 
for an initial water temperature of 20 °C is the maximum value for which evaporation is stil] 
possible and beyond which nucleate boiling is observed. By setting the temperature ratio 
in the figure to a common value (= 0.8), similar interfacial conditions are expected for the 


various materials. Recall the temperature profiles of Fig. Al which clearly indicate the 
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different local temperature drop for aluminum and Macor. This explains why similar local 
temperatures at the droplet outer edge are observed for bakelite with an initial surface 
temperature of 230 °C and for aluminum with an initial surface temperature of 88 °C. Note 
the temporal behavior of the spatial averaged heat flux which in all cases tends to sharply 
rise toward the end of the process due to the thinning of the water layer and to the 
subsequent sharp increase in heat transfer. 

Figure A8 shows the corresponding temperature profiles toward the end of the evaporation 
transient (t/r = 0.9). Note the dip in the temperature near the edge of the droplet for the 
paeliee which, due to its very low conductivity, enhances the local temperature gradients. 
In all cases note that the temperature at the edge of the droplet is higher than the spatial 
averaged value. Also note the different behavior of the surface temperature distributions 
outside the wetted area and in particular the radial distance where the effect of the droplet 
cooling becomes negligible. 

Solid-liquid interfacial temperature profiles for Macor and aluminum are shown in Fig. A9. 
The high thermal conductivity of aluminum results in uniform and constant interfacial 
temperatures while the lower thermal conductivity of Macor is reflected in non-uniform and 
variable temperatures. The Macor interfacial temperatures confirm the discrepancies 
observed in Fig. A5 and demonstrate the limits of the assumptions used by diMarzo and 
Evans (1989), by Seki et al. (1978) and by Rizza (1981) in their simplified evaporative 


cooling analyses. 
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A4.3. Evaporation time correlation 

One of the most significant parameters in the analysis of the results is the total evaporation 
time. The lack of approximate solutions requires the formulation of a correlation. A simple 
form for the normalized evaporation time is fitted with the experimental data for aluminum 


and Macor yielding the following result: 


po Lee (aazy° (=| x 


R B a, 
(A17,A18) 
k.\0-640 7 kc. \9-0858 
e -2.45]—= See Fee. |e= 
| (=) T,-T, ; | 


The comparison between the total evaporation times provided by this correlation and those 
predicted by the model for various materials and for various conditions is shown in Fig. A10. 
The correlation is in excellent agreement with the model for evaporation times below 100 
seconds. For higher evaporation time it tends to over-predict within 20 per cent of the 
model results. It is notable that the correlation is able to reasonably describe the 


evaporation time with the normalized parameters identified by the approximate solution. 


A4.4. Radius of influence 

In order to quantify the extent of the solid surface region affected by the droplet cooling, 
a significant parameter is the radial heat flux at the surface at any given location. The 
presence of the droplet is felt when this radial heat flux is greater than a reference value 
based on the spatial and temporal averaged heat flux, that is: 
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hope. eee Ap,V 
Ke ar 2 ty) aia (A19) 


By using Eq. (A13) and Eq. (A14), this condition becomes in terms of normalized quantities: 


a. 4 
ak ae 


an 
fo AGn) a(S) ext (68) dz 


(A20,A21) 


where ¢ and yf are arbitrary constants. Equations (A20) and (A21) define the radius of 
influence n, (= 1,/R) as the minimum value of 7 that satisfies these inequalities. Two 
different constants are used because the approximate closed-form solution, which yields Eq. 
(A21), is able to capture the trends of the radius of influence but not the actual value as 
clearly seen in Fig. All. The radius of influence can be defined for @ = 0.1 which relates 
the radial heat flux to ten per cent of the spatial and temporal averaged heat flux [see Eq. 
(A19)]. The approximate, closed-form solution for y = ¢ = 0.1 is shown in the figure. 
wales of y between 0.045 and 0.065 provide a good representation of the model predictions 


for a wide variety of conditions. 


The combination of Eq. (A18) and of Fig. All provides a rapid estimate of the radius of 
influence. Both the evaporation time correlation and the radius of influence plots are 
excellent predictive tools at low values of the normalized time 6. The usual values of this 


parameter for most materials and conditions is found to be below five. 
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A5. CONCLUSIONS 


This study briefly reviews the formulation of a model for the prediction of the thermal 
behavior of the droplet-solid interaction during evaporative cooling. The model is validated 
over a wide range of parameters. Model predictions for a variety of conditions are used to 


gain insight into the evaporative process mechanism. 


A correlation for the prediction of the total droplet evaporation time is proposed and 
checked against the model predictions. An approximate, closed-form solution is used to 
normalize the relevant parameters and to provide estimates of a newly defined radius of 


influence. 
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Aé. 


Zepeda alee ae oe 


i 
7) 


NOMENCLATURE 


thickness of the droplet 
specific heat 


steam water mass diffusivity - diameter 


error function 

forcing and unknown vector 
Green’s function 

overall heat transfer coefficient 
convective heat transfer coefficient 
latent heat of vaporization 

Bessel’s functions 

Jakob number 

thermal conductivity 

modified Bessel’s function: e* I,(7) 
Lewis number: (D/a)** 

unit vector norma! to the surface 
reference heat flux: A p V / (x R’r) 
spatial averaged heat flux 

radial coordinate ~ 

radius of the wetted area 

surface 

time 

recollection time 

temperature 

onset of nucleate boiling temperature 
contact temperature 

transformed temperature 

volume or drop volume 

weight matrix 

molar fraction of steam in air 


axial coordinate 
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thermal diffusivity 

shape parameter 

normalized time: (a, t)/?/R 

solid surface emissivity 

dummy variables 

normalized radius: r/R 

radius of influence associated with the constant ¢ 
normalized temperature 

liquid latent heat of vaporization 


‘density 


Stefan-Boltzmann constant 
total evaporation time 
arbitrary constants 


ambient, air 

liquid-vapor interfacial property 
liquid i 
generic point index 

in the radial direction 

solid 

in the axial direction 

Green’s function argument 
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Fig. Al Typical solid surface temperature profiles for water droplets (V = 30 ul, r = 
95 s) deposited on aluminum (T, = 91 °C) and Macor (T, = 119 °C) for a) 
t/r = 03 and b) t/r = 0.9 7 : 
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Fig. A2 
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Fig. A3 Model validation: solid surface temperatures for water on Macor at t/t = 0.9 


with a) V = 30 ul, T, = 160°C; b) V = 30 wl, T, = 101°C; c) V = 10 ul, T, 
= 143 °C and d) V = 30 wl, T, = 143 °C (data from Klassen et al. 1990) 
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Fig. A4 
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Model validation: total evaporation time for droplets with V = 10, 30, 50 ul 
on aluminum with T, = 75 + 102 °C and on Macor with T, = 101 + 208 °C 
(data from Klassen et al. 1990) . 
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Approximate, closed-form solution versus experimental data for water droplet 
on Macor with a) V = 30 ul, T, = 160°C; b) V = 30 wl, T, = 101°C; c) V 


= 10 ul, T, = 


1990) 


143 °C and d) V = 50 ul, T, = 143 °C (data by Klassen et al., 
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Fig. A6 
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Liquid-solid interfacial fluxes for water droplet with V = 30 ul on Macau with 
Tess 143 °C at d) t/r = 0.3; c) t/r = 0.5; b) t/r = 0.7; a) t/r = 0.9 and r 
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Normalized spatial averaged liquid-solid interfacial flux for a water droplet 
with V = 30 ul on a) bakelite with T, = 230 °C; b) aluminum with T, = 88 
°C; c) marble with T, = 125 °C; d) steel with T, = 98 °C; T, from Eq. (1) and 
T, = 20°C 
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Typical solid surface temperature profiles for a water droplet with V = 30 ul 
at t/t = 0.9 on a) bakelite with T, = 230 °C; b) macor with T, = 143 °C; c) 
marble with T, = 125 °C; d) steel with T, = 98 °C; e) aluminum with T, = 88 
°C; T, from Eq. (1) and T, = 20°C 
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Fig. A9 Typical liquid-solid interfacial temperatures for Macor and aluminum with T, 
= 82°C for water droplets with V = 30 ul at a) t/r = 03 and b) t/r = 0.9 
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Fig. A10 
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Evaporation time correlation versus model predictions for water droplets with 
V = 10,-30, 50 ul on various materials with T, = 80 + 230 °C 
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Radius of influence versus normalized time: closed-form, approximate solution 
and model predictions where ¢ = y = 0.1; a) ¢ = 0.065 and b) # = 0.045 


Table 1 - PROPERTIES OF ALUMINUM AND MACOR 


ALUMINUM 


Density: 


Thermal Conductivity: 


Specific Heat: 


RESEKLELLSE 


MACOR 


Density: 


Thermal Conductivity: 


Specific Heat: 


2771 kg/m? 
180 W/(mK) 


962 J/(kgK) 


2520 kg/m? 
1297 W/(mXK) 


888.9 J/(keK) 
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PART B 


RADIATION EXPERIMENTS 


Bl. INTRODUCTION 


Quenching of metals, spray cooling and fire suppression with sprinklers are few of the 
applications that inspired numerous theoretical and experimental investigations. Global 
approach to spray-surface interactions is addressed in the early investigations by Toda [1] 
and Bonacina [2]. More detailed numerical and experimental studies followed. In 
particular, the contributions of Rizza [3], Seki et al [4], Inada et al. [5], diMarzo and Evans 
[6][7] are based on the integration of the transient conduction equation in the solid with a 
simple boundary condition at the liquid solid interface which prescribes constant and 
uniform temperature throughout the vaporization process. Grissom: and Wierand [8] 
determined experimentally the lowest surface temperature possible for the existence of spray. 
evaporative cooling and presented a conduction-controlled analytical model of droplet 


evaporation. 


Other experimental works reinforced these modeling assumptions for solids characterized 
by high thermal conductivity. In particular, Sadhal and Plesset [9] described the droplet 
shape behavior when the liquid thins and breaks up; Pedersen [10] presented a large 
number of experimental data, showing that, for impinging droplets, approach velocity is a 
dominant variable affecting droplet heat transfer; Choi and Yao [11] studied the heat 
transfer characteristics of impacting sprays, with special attention to the differences between 


sprays impacting horizontally and vertically. Makino and Michiyosi [12][13][14][15] 


B-2 


investigated the onset of nucleate boiling and the corresponding behavior of the droplets 
during the different stages of the evaporating process, with particular interest on the 
influence of the governing parameters such as the droplet size and the initial surface 
temperature. These investigations span the full range of vaporization phenomena to include 
evaporation, nucleate boiling and film boiling. Detailed information on the droplet shape 


[9][13] and on the solid thermal transient [12][14][15] are given. 


Recent studies by Abu-Zaid and Atreya [16], diMarzo and Evans [6][7] and Klassen et al. 
[17] concentrated on the evaporative phenomena with temperature ranges below the onset 
of nucleate boiling. Low thermal conductivity solids are heated from below and the 
evaporation of a single droplet is investigated [6][7][17]. The behavior of porous and non- 
porous materials with similar thermal properties are compared [16] and porous materials 


show higher evaporative cooling rates. 


This work is part of a long term research effort aimed at modeling the extinguishment of 
a solid fuel fire. Water, used as an extinguishing agent, is a primary cause of secondary fire 
damage. By understanding and predicting fire extinguishment mechanisms, one can optimize 
the water application (droplet sizes, mass fluxes) and thus minimize the damages associated 
with the extinguishment process. Previous studies focused on the evaporation of water 
droplets deposited on surfaces which were heated by conduction from below. Klassen [17] 
studied droplet evaporation from Macor (a glass-like material) and diMarzo [6][7] reported 


on the evaporative cooling of aluminum. Most of these data refer to the range where only 
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evaporation takes place and nucleate boiling is fully suppressed. Some data points are also 
taken in the nucleate boiling range. In the present study, the heat input is provided by 
radiant panels located above the Macor to simulate a fire environment. Data are compared 
with Klassen’s previous experiments [17] to determine the differences in the evaporation 
process due to the heat input characteristics (i.e. radiation from above versus conduction 
from below). Klassen [17] developed an infrared thermographic technique to obtain 
transient temperature profiles on the solid surface during the droplet evaporation. This 
technique is non-intrusive and has no response lag. The primary purposes of this work are: 
a) to gather data which will be used to validate models of the evaporation phenomena and 


b) to provide insight into dropwise evaporative cooling. 
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B2. EXPERIMENTAL APPARATUS 


The experimental apparatus (shown schematically in Fig. 1) consists of a Macor tile (of 
squared shape with 15.2 cm sides and 2.54 thick) mounted onto a steel chilled plate. 
Contact between the tile and the chilled plate is enhanced by a conductive paste. The 
Macor is located under two radiant panels facing each other and oriented towards the 
Macor surface. An infrared camera is mounted above the tile at a 45 degree angle from the 
vertical axis and looks through a chilled pipe at the Macor surface. A mechanized droplet 
dispenser rolls in and out of the heated region and is designed to dispense metered water 


volumes at a precise location and in a repeatable manner (see for details Ref. [6]). 


The radiant panels are used to simulate the heat flux output of a fire. These panels are 
powered by a three-phase 220 volt power supply which is controlled by a 0-10 volts output 
temperature controller. The panel surface temperatures provide the feedback to the 


controller. 


A chill plate is placed under the Macor tile and is kept at near ambient temperature by 
controlling the water flowing through it. The heat flux in the vertical direction overwhelms 
the heat losses through the side of the tile. This set up is designed to obtain (prior to the 
droplet deposition) a uniform Macor surface temperature by reducing the radial 


temperature gradients and by providing a linear temperature profile across the tile in the 
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vertical direction. 


The infrared camera (Model 525 Infrared Thermal Imaging System by Inframetrics) is used 
to observe the temperature distribution on the Macor tile surface. The camera is connected 
to a VCR which records the transient infrared temperature distribution of the portion of the 
Macor surface affected by the water droplet. The video image is recorded in real time. A 
chilled pipe is used to minimize the amount of extraneous infrared radiation that reaches 
the camera lens. This chilled pipe is constructed with a water cooled copper coil. The 
interior of the chilled pipe is coated with a high emissivity black paint to minimize the 


amount of reflection. 


A mechanized droplet dispenser is used to place the water droplets onto the Macor surface. 
The dispenser rolls in and out of the experimental apparatus quickly and easily so that 
almost no disruption of the radiant heat flux is observed. A microdispenser (Drummond 
Scientific Corporation 550) is used to measure the droplet volume and to dispense it on the 
Macor surface. The dispenser is operated by a motor driven cam to obtain a repeatable 


droplet shape on the surface after deposition. 


Macor properties are listed in Table 1. In particular note that it has high emissivity, which 
minimizes reflection from the surface. Macor can withstand repeated thermal shocks and 


its surface is smooth and cracks free. 


B3. EXPERIMENTAL PROCEDURES 


An extensive test matrix is implemented following an identical experimental protocol. These 
experiments are conducted for droplet sizes of 10, 30 and 50 ul at surface temperatures of 
90, 110, 130, 145, 150, 155,.and 164 °C. The droplet sizes are a compromise between the 
sprinkler delivered droplets (= 1 yl) and reasonably large volumes to achieve easy 
repeatable dispensing (> 10 ul). Note that surface tension impairs the delivery of smaller 
single droplets under the sole gravitational forces. The temperatures are mostly in the 


evaporative range since onset of nucleate boiling is observed at and above 155 °C. 


The experimental protocol consists of 5 steps: 


B3.1. Water conditioning 

The water used in the experiment is de-ionized and must be completely degassed. De- 
ionized water is used to maintain a clean, deposit-free Macor surface. Complete degassing 
must be achieved to avoid degassing during the evaporative process. If gasses come out of 
solution during the evaporation they may alter the droplet shape and therefore its behavior 


or they can even shatter the droplet or trigger nucleate boiling. 


Degassing is obtained by repeatedly freezing and thawing water kept under vacuum. The 


water is then removed from the degassing container by bringing it to boil and by extracting 
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the liquid with a syringe. The syringe is then used to load the mechanized dispenser. 


B3.2. Infrared camera, radiant panels, chilled plate and pipe 

The radiant panels and the chilled plate are activated and the Macor tile is brought to 
thermal steady state. The infrared camera sensing enclosure is constantly monitored to 
insure that adequate cooling by liquid nitrogen is available (a liquid nitrogen charge lasts 
about 45 minutes). The chilled pipe is also activated and the infrared picture of the Macor 
surface is obtained. The infrared camera provides a two-dimensional picture where the 
vertical axis relates to a radiance reading while the horizontal axis is the actual position of 
the radiance reading on a line which can be arbitrarily set on the Macor surface. This line 
is carefully positioned through the center of the wetted area (i.e. the liquid-solid interfacial 


area). 


B3.3. Temperature calibration 

The radiant panels control set point is adjusted and adequate time is allowed for the system 
to reach thermal steady state. The infrared picture needs to be calibrated to relate the 
radiance readings to temperatures. The temperature scale is provided by the infrared 
camera manufacturer and is also calibrated independently by using a blackbody source. The 
temperature scale is referenced with one radiance/temperature point. This reference point 
is obtained by measuring the surface temperature with a hand held probe (QMEGA surface 
probe series 68000) until a constant reading is observed. By quickly removing the probe 


while recording the infrared picture of the probe location, one obtains an instantaneous 
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reading of the same location with the infrared camera and with the temperature probe. 
Note that the temperature at the probe location changes after the probe removal and that 
only the picture immediately after the probe removal provides the desired information. 
Since 30 frames per second are acquired with the VCR, there is no problem in selecting the 


appropriate frame. 


B3.4. Droplet deposition 

Prior to the droplet deposition, the surface is cleaned with alcohol and with de-ionized 
water. When the cleaning water is evaporated the prope! is deposited. The line on which 
the temperatures are scanned is adjusted to pass exactly through the center of the wetted 
region. This is a fine tuning adjustment since the mechanized dispenser places the droplet 
within fractions of a millimeter of the desired location. The evaporative process is then 
recorded and set aside for post-processing. Each experiment is repeated at least ten times 


to build sufficient confidence in its repeatability. 


B3.5. Test record 

A regular black and white camera is used to record (on the same VCR tape) the 
experimental conditions. In particular, droplet size, Macor surface temperature, chilled 
plate temperature and radiant panels temperatures are recorded. In addition, an infrared 
picture of a small coin placed on the Macor surface is also recorded to determine the length 


scale of the infrared image. All this information is recorded at the beginning and at the end 
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of a ten tests sequence (which relate to the same droplet size and initial temperature 


conditions). 
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B4. DATA PROCESSING 


After a sequence of ten tests is completed, data processing is performed. The evaporation 
time + for each experiment is measured from the VCR record. The size of the wetted 
region is identified and the droplet shape is inferred. Extensive visual observations and 
photographic records show that the droplet shape can be approximated by a segment of a 
sphere [6][9]. By inspecting the droplet shape and evaporation time information, one 
representative test is selected out of the ten tests sequence. The rationale for this selection 


is based on the following concepts: 


a) The droplet evaporation time increases as solid-liquid interfacial area decreases. 
Since the droplet geometry is reasonably approximated by a segment of a sphere, a 
shape parameter, 8 [2][6][9], is defined as the radius of the wetted area over the 
radius of an equal volume droplet of spherical shape. An averaged and most 


repeatable shape parameter identifies the representative droplet. 


b) In spite of the dispenser accuracy, some of the water will cling to the inside of the 


droplet dispenser. The longer a droplet takes to evaporate, the greater the likelihood 


te 


that less water is left in the dispenser. 


c) The droplet does not always fall in the same location. Sometimes the temperature 
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scanned line is not positioned exactly through the droplet center. Further, the wetted 
region shrinks during the evaporation and the center of the wetted region may drift 
from its initial position. If this occurs, the image may show complete evaporation 


when water is still present. 


These three concepts are the guidelines for the selection of the most representative droplet. 
Usually there are four or five identical tests. In a few cases, a decision must be made based 


on item a) with some bias towards higher evaporation times to account for items b) and c). 


Images from the record for the selected droplet are obtained at 10, 30, 50 and 100 per cent 
of the evaporation time and a few images are also analyzed to describe the surface 
temperature recovery after the evaporation. These images are digitized by a PC-based 641H 
x 201V frame grabber (by Computer Eyes which also provides the required software). The 
spatial resolution in the horizontal axis may vary depending on the camera field of view. 
Usually a pixel corresponds to 0.05 millimeters. The temperature scale is of about 100 °C 
and the typical thermographic trace is two pixels thick. Therefore, a resolution of 1-2 °C 


is achieved in the vertical axis (i.e. the temperature scale). 


The digitized image is processed with a thresholding and subsequent erosion techniques 
which are described at length by Klassen [17]. In summary, the intensity of the 
thermographic traces is differentiated from unwanted noise by thresholding. This technique 


identifies as noise the pixels with intensities above a set value and retains the pixels with 
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intensity below this set threshold. In addition, thresholding binarizes the resulting image by 
assigning the maximum and minimum intensity value respectively to data pixels or "blank" 
pixels. Further elimination of noise is achieved by erosion which consists in a pixel-by-pixel 
inspection of the neighboring information to differentiate between isolated information 
(which is most likely noise) and clustered information (which is most likely useful data). 
These techniques are applied by trial and error and visual inspection of the images confirms 


their effectiveness. 


The data points isolated by the image processing techniques are then compacted by folding 
the right side of the image (about the droplet center) onto the left side of the image. Since 
the temperature profile is symmetric about the droplet center, superimposing one half of the 
profile on the other half produces twice as many data points to be analyzed and allows a 
more accurate determination of the symmetry axis through the center of the wetted region. 
Figure 2 shows the typical transient data after processing. Note that in the figures prior to 
the complete evaporation (t < 1), the temperature profiles are plotted for r 2 R. The 
radiance signal in the solid surface region covered by the droplet cannot be translated into 
a significant temperature reading since water is a radiation participating medium. The 
temperature distributions in the range 0.3 r < t s 0.9 rt do not change appreciably and the 


surface temperature recovers its initial value for t 2 1.5 Tr. 


B5. RESULTS AND DISCUSSION 


B5.1. Phenomenology 

To set the results and their discussion in the appropriate context it is necessary to briefly 
review the major differences in the thermal behavior between the case of radiant heat input 
(from above the surface) and the case of conductive heat input (from below the surface). 
It is important to note that, for the droplet sizes under investigation, the heat transfer 
mechanisms in the liquid are conduction and radiation. Convection does not play a 
significant role since Rayleigh-Benard instabilities require times longer than the evaporative 
time to develop and Marangoni surface tension driven flows generate very slow motion in 


these experimental conditions (Ostrach and Pradhan [18]). 


In the case of radiant heat input, the initial vertical temperature gradient in the solid shows 
a higher temperature at the surface and a lower temperature at the chilled plate. The 
droplet tends to reduce the heat transfer in the vertical direction since it decreases the 
vertical temperature gradient. In the case of conductive heat input, the initial vertical 
temperature gradient is opposite since the temperature increases with depth. Droplet 
cooling increases the heat transfer in the vertical direction since it increases the vertical 
temperature gradient. This means that the portion of the solid under the droplet 
contributes greatly to the evaporation of a droplet in the conduction case while it has a 


lesser role in the radiation case. Note that for both cases the amount of heat required to 
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evaporate a given size droplet is the same. 


The conductively heated surface is certainly supplying more energy to the droplet through 
conduction than the radiatively heated surface, yet Fig. 3 shows that the radiatively heated 
surface evaporation times are shorter. This is due to the heat directly supplied by radiation 
to the water. The radiant panels behave like blackbodies. Their typical surface temperature 
is in the 600-750 °C range. At these temperatures the maximum emission of the blackbody 
spectrum occurs in the 3-3.5 micron range (Siegel and Howell [19]). Water has its maximum 
absorption coefficient of 11,400 cm”? at a wavelength of 3 microns and maintains an elevated 
value at all wavelengths in the near infrared above 2.6 microns (Siegel and Howell [19]). 
The absorption coefficient determines the exponential attenuation of the light intensity. By 
multiplying the absorption coefficient by the thickness of the water layer, one obtains a non- 
dimensional exponent which determines how much radiant energy is transmitted through the 
water. Figure 4 shows the percentage of incident radiation which is absorbed by water (10, 
30 and 50 yl droplets) at various initial solid surface temperatures, corresponding to various 
radiant panels temperatures. The line in the figure is the linear regression of the data. The 
temperature distribution at the exposed solid surface is shown in Fig. 5 for both the 
radiation and the conduction case. The curves a) and b) represent respectively a radiation 
and a conduction case with the same initial surface temperature, while curve c) shows the 
temperature distribution for a radiation case in which the evaporation time is the same of 
the conduction case b). All curves are obtained for a 30 ul droplet. These three curves are 


for t = 0.3 1 since thereafter the droplet shrinks in the radiation case and a graphical 
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representation of the temperature profiles would be quite intricate due to the time 


dependence of the normalizing parameter R. 


A direct comparison of the radial heat flux is not possible due to the uncertainty associated 
with the data at the droplet edge (r = R). However it is apparent that the cooling effect 
is less pronounced for the radiation case since radiation contributes directly to the 


vaporization process thus reducing the conduction component. 


BS5.2. Shape parameter 

Figure 6 shows the shape parameter, £, for various droplet sizes and initial solid surface 
temperatures. The shape parameter is defined as 8 = R/R,; where R is the radius of the 
wetted region under the deposited droplet and R, is the radius of a spherical droplet of the 
same volume. Note that an experimental measurement of f is possible only at the 
beginning of the evaporative transient, because it is not possible to obtain information about 
the droplet volume during the evaporation, thus R, cannot be determined after the initial 


time t=0. 


The droplet size is not significantly affecting 8. This holds true for the conduction case. 
The scatter is overwhelming for temperatures above 155 °C. This temperature corresponds 
to the onset of nucleate boiling. For the conduction case, nucleate boiling is observed at 


164 °C and above. Klassen [17] measured @ values that ranged from 1.2 to 1.5. These 
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values are significantly lower than the values measured for the radiation case. Klassen also 
showed constant radii of the wetted region. Similar results were obtained by diMarzo [6] 
for a quartz surface. Both these results were obtained for the conduction case. It appears 
that, for the conduction case, increasing the solid stds temperature does not affect the 
water-vapor interfacial temperature. This is reasonable since the heat can only be 


conducted from the solid surface through the liquid layer. 


For the radiation case, B average values range from 1.4 at an initial solid surface 
temperature of 90 °C to 2.5 at an initial solid surface temperature of 155 °C. This may be 
due to the direct radiation from above the droplet. The absorption coefficients are high for 
water in the spectra produced by the radiant panels. This results in substantial amounts of 
energy absorbed by the surface of the water, causing the droplet surface to heat up and 
therefore reducing its surface tension. This in turns may cause the droplet to spread out. 
The higher is the temperature of the heaters, the greater is the amount of heat absorbed 
by the droplets, which lowers the surface tension and contributes to a wider spreading of the 


droplet. 


Figure 7 shows how the radius of the wetted region changes during the evaporation process 
for various initial solid surface temperatures. Note that the radius of the wetted region R 
is normalized with its value at deposition R,. At higher initial surface temperatures (i.e. at 
the onset of boiling), a marked decrease in the radius of the wetted region is observed as 


the evaporation progresses. This shrinking of the wetted region occurs when the solid-liquid 
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contact angle has reached its receding value (Simon and Hsu [20]). This phenomenon is 


observed by Klassen [17] at the very end of the evaporative transient. 


Although there are many differences between the solid thermal response for the conduction 
and radiation cases, there are also some similarities. Both cases show an onset of nucleate 
boiling when the initial surface temperature is in the 155-164 °C range. Both cases show 


that nucleate boiling drastically disrupts the droplet geometry with large variations of £. 
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B6. CONCLUSIONS 


The experiments describe some aspects of the behavior of a low thermal conductivity 
material heated by radiation from above when a water droplet is placed on its surface. The 
results are different from analogous results obtained with the same material heated by 


conduction from below. 


Water absorbs near infrared radiation. More than 95 per cent of the radiant heat supplied 
from the radiant panels is in the near infrared spectrum. It is estimated that 20 to 30 per 
cent of the heat of vaporization is supplied to the water droplet from direct absorption of 


the incident radiation. The remainder is from conduction across the liquid-solid interface. 


The shape parameter increases as the surface temperature increases. A contributing factor 
is the effect of surface tension at different temperatures. For the conduction case, the 
constant value of @ is an indicator that the temperature at the exposed water droplet surface 


is not affected by the initial solid surface temperature. 


The evaporation times for the radiation case are shorter than for the conduction case. The 
water droplets on the radiatively heated surface are heated from both below and above. 
Further, the radiatively heated droplets spread across the surface more than in the 
conduction case thus providing a larger surface area to absorb heat from above and a 


thinner layer to transfer heat from below by conduction. 
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ot 


NOMENCLATURE 

radial coordinate on the solid surface 

radius of the wetted region [f(t) in the radiation case] 

radius of the wetted region at deposition 

radius of a sphere with equal volume of the droplet at deposition 
time 

initial solid surface temperature 


volume of the droplet at deposition 


shape parameter: = R/R, 


total evaporation time 
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FIGURE 1 - Experimental apparatus (top view) 
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FIGURE 2 - Transient temperature distribution (T, = 130°C; V = 30 ul; 7 = 605) 
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FIGURE 3 - Evaporation time (solid symbols for radiation case; open symbols for 


conduction case; o: V = 50 ul; a: V = 30 wl; O: V = 10 pl) 
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FIGURE 4 - Percentage of vaporization heat by direct radiation ( a: V = 50 ul; o: V = 30 


ul; O: V = 10 pl) 
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FIGURE 5 - Cooling effect: comparison between radiation and conduction (V= 30 ul; a: 
radiation T, = 145 °C; r = 40s; D: conduction Ty = 145 °C; r = 605; c: 


radiation T, = 130°C; r = 60s ) 
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SHAPE FACTOR 
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FIGURE .6 - Shape factor (o: V = 50 wl; a: V = 30 wl; O: V = 10 yl) 
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FIGURE 7 - 


~ INITIAL SURFACE TEMPERATURE (°C) 


Normalized radius of the wetted region (0: t = 0.3 3; a: t = 0.6 7; 0: t = 0.9 


T) 


ASR REET SEAT ed EI IS MEET PE BSE SEE TI 
TABLE 1 - MACOR PROPERTIES 


Density: 2520 kg/m? 
Thermal Conductivity: 1.297 W/ (m°K) 
Specific Heat: © 888.9 J/ (kg°K) 
Emissivity: 0.94 


B- 32 


B10. APPENDIX - DATA BASE 


In the following pages, experimental data from all the tests conducted are presented. 
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INITIAL INIT. TEMP. CHILL PLATE EVAPORATION BETA 


VOL. (pl) MACOR (°C) TEMP. (°C) TIME (s) (adim.) 
10 180 59 8.41 3.39 
10 180 59 725 3.81 
10 180 59 a7) 3.47 
10 180 59 6.46 3.96 
10 180 59 6.98 3.88 
10 180 59 736 3.12 
10 180 59 5.81 3.24 
10 180 59 7.03 3.64 
10 180 59 6.59 3.58 
10 180 59 6.59 3.29 
30 180 59 10.91 3.86 
30 180 59 11.16 3.34 
30 180 59 3 9.65 3.34 
30 180 59 9.53 3.36 
30 180 59 10.46 3.11 
30 180 59 8.64 3.27 
30 180 59 8.16 3.29 
30 180 59 10.47 3.32 
30 180 59 9.34 3.84 
30 180 59 9.15 3.92 
30 180 59 12.42 3.15 
50 180 59 | 12.42 3.34 
50 180 59 14.30 3.26 
50 180 59 8.05 3.57 
50 180 59 11.81 3.84 
50 180 59 13.29 3.06 
50 180 59 10.97 3.19 
50 180 59 14.36 3.18 
50 180 59 5.56 239 
50 180 59 12.62 3.15 
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INITIAL INIT. TEMP. CHILL PLATE EVAPORATION BETA 


VOL. (pl) MACOR (°C) TEMP. (°C) TIME (s) (adimn.) 
10 164 37 11.90 3.06 
10 164 37 13.59 2.98 
10 loa’ 37 15.30 2.94 
10 164 37 11.84 1321 
10 164 37 12.86 3.39 
10 164 37 15.53 2.88 
10 164 37 5.79 2.41 
10 164 37 13.00 2.86 
10 164 37 , 13.90 3.04 
10 164 37 13.09 3.09 

30 164 37 24.72 2.60 
30 164 37 17.84 4.02 
30 164 37 15.22 4.12 
30 164 37 23.79 3.03 
30 164 37 23.98 3.03 
30 164 37 30.99 2.72 
30 164 oy 16.59 3.23 
30 164 37 | 27.03 2.81 
30 164 37 29.55 pep 
50 164 37 22.15 2.98 
50 164 37 29.84 3.04 
50 164 37 32.00 3.29 
50 164 37 29.98 3.41 
50 164 37 20.28 2.98 
50 164 37 33.09 2.76 
50 164 37 28.73 3.21 
50 164 37 22.34 3.47 
50 164 37 20.36 3.19 


50 164 nw 25.47 Path 
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INITIAL INIT. TEMP. CHILL PLATE EVAPORATION BETA 


VOL. (jl) MACOR (°C) TEMP. (°C) TIME (s) (adirn.) 
10 155 36 14.87 2.88 
10 155 36 16.03 2.59 
10 me 155 36 16.86 3.02 
10 155 36 14.15 3.45 
10 155 36 16.96 | 2.98 
10 155 36 16.53 3.08 
10 155 36 14.75 3.02 
10 155 36 12.59 3.14 
10 155 36 1653 3.08 
10 155 36 16.59 2.86 
30 155 36 28.44 2.87 
30 155 36 33.65 2.60 
30 155 36 33.59 2.80 
30 155 36 3121 2.64 
30 155 36 27.96 2.93 
30 155 36 32.21 2.54 
30 155 36 31.30 2.57 
30 155 36 32.37 2.60 
30 155 36 32.85 2.52 
30 155 36 32.80 2.57 
50 155 36 39.90 2.47 
50 155 36 43.00 2.38 
50 155 36 40.90 2.44 
50 155 36 38.03 2.42 
50 155 36 : 41.03 2.73 
50 155 36 38.55 2.58 
50 155 36 37.22 230 
50 | 155 36 42.93 2.36 
50 ee) alg 36 45.97 2.41 
50 155 36 45.28 2.39 
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INITIAL INIT. TEMP. CHILL PLATE EVAPORATION BETA 


VOL. (ul) | MACOR (°C) ‘TEMP. (°C) TIME (s) (adirm.) 
10 150 35 16.46 2.55 
10 150 35 20.15 2.60 
10 150 35 14.97 2.43 
10 150 35 14.70 2.46 
10 150 35 19.58 2.55 
10 150 35 17.42 2.62 
10 150 35 14,05 2.61 
10 150 35 17.89 2.55 
10 150 35 18.25 2.44 
10 150 35 16.28 2.49 
30 150 35 32.90 2.45 
30 150 35 32.44 | 2.32 
30 150 35 29.85 2.54 
30 150 35 28.46 2.56 
30 150 35 33.46 2.40 
30 150 35 35.42 2.38 
30 150 35 28.47 2.50 
30 150 35 34.53 2.30 
30 150 35 31.40 2.41 
30 150 35 32.24 - 240 
50 150 35 38.74 2.47 
50 150 35 39.36 2.40 
50 150 35 42.53 ese 
50 150 35 38.36 254 
50 150 35 37.30 2.32 
50 150 35 39.22 2.38 
50 150 35 33.85 2.46 
50 150 35 38.66 251 
50 150 35 39.09 2.30 
50 150 35 36.41 2.33 
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INITIAL INTT. TEMP. __ CHILL PLATE 


VOL. (ul) 


10 
10 
10 
10 
10 
10 
10 
10 
10 
10 


30 
30 
30 
30 
30 
30 
30 
30 
30 
30 


50 
50 
50 
50 
50 
50 
50 
50 
50 
50 


MACOR (°C) 


145 
145 
145 
145 
145 
145 
145 
145 
145 
145 


145 
145 
145 
145 
145 
145 
145 
145 
145 
145 


145 
145 
145 
145 
145 
145 
145 
145 
145 
145 


TEMP. (°C) 


35 
Sh) 
35 
35 
35 
35 
35 
33 
35 
35 


35 
o2 
55 
55 
35 
23 
35 
35 
ape) 
35 


35 
35 
35 
35 
35 
Sh 
a0 
35 
35 
sh) 


EVAPORATION 
TIME (s) 


21.78 
10°50 
22.86 
17.34 
20.87 
21.47 
19.34 
20.03 
19.42 
21.54 


38.36 
40.91 
40.40 
36.91 
37.06 
40.65 
35.96 
39.06 
38.37 
41.35 


48.10 
49.09 
33.39) 
46.00 
52.80 
51.03 
47.06 
50.28 
51.92 
32.97, 


INITIAL INIT. TEMP. CHILL PLATE EVAPORATION BETA 


VOL. (pl) MACOR (°C) TEMP. (°C) TIME (s) (adim.) 
10 130 34 30.40 1.96 
10 230 34 31.81 1.93 
10 130 34 32.00 2.01 
10 130 34 32.15 | 2.05 
10 130 34 30.59 ek 93 
10 130 34 33.15 1.91 
10 130 34 29.71 1.89 
10 130 34 29.30 1.93 
10 130 34 31.12 2.01 
10 130 34 31.90 2.08 
30 130 34 55.46 1.96 
30 130 34 59.71 1.94 
30 130 34 59.96 1.93 
30 130 34 54.40 1.88 
30 130 34 55.87 1.84 
30 130 34 5721 1.90 
30 130 34 54.03 1.94 
30 130 34 58.90 1.93 
30 _ 130 34 58.96 1.90 
30 130 34 59.15 1.93 
50 130 34 87.00 1.79 
50 } 2150 34 83.71 1.78 
50 130 34. 86.21 1.82 
50 130 34 84.04 1.84 
50 130 34 81.53 1.91 
50 130 34 84.43 1.83 
50 130 34 84.53 «1.84 
50 130 34 85.97 1.90 
50 130 34 84.31 1.80 
50 130 34 85.35 1.83 
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INITIAL INIT. TEMP. CHILETIATE . EVAPORATION BETA 


VOL. (ul) MACOR (°C) TEMP. (°C) TIME (s) (adim. ) 
10 110 33 45.21 1.81 
10 110 33 | 45.90 1.83 
10 eit 33 46.28 1.75 
10 110 33 46.59 1.73 
10 110 33 47.42 1.80 
10 110 33 46.71 1.79 
10 110 33 45.53 1.75 
10 110 33 44.84 1.82 
10 110 33 46.34 1.73 
10 110 33 46.71 1.80 
30 110 33 83.60 1.69 
30 110 33 85.00 1.73 
30 110 33 84.23 1.76 
30 110 33 81.23 | 1.81 
30 110 33 82.16. 1.70 
30 BID so 33 85.67 1.71 
30 110 33 72.53 1.68 
30 110 33 84.00 1.70 
30 110 33 84.24 1.71 
30 110 33 83.40 1.70 
50 110 33 110.21 1.73 
50 110 33 110.84 1.71 
50 110 33 110.92 1.69 
50 110 33 112.60 1.71 
50 110 33 110.21 1.65 
50 110 33 110.97 1.69 
50 110 33 107.84 1.70 
50 moO. 33 101.96 1.71 
50 110 33 102.59 1.68 
50 110 33 109.53 1.70 
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INITIAL INIT. TEMP. CHILL PLATE EVAPORATION BETA 


VOL. (ul) MACOR (°C) TEMP. (°C) TIME (s) (adim.) 
10 90 30 75.66 1.63 
10 90 30 75.40 1.62 
10 . 90 30 76.30 1.62 
10 90 30 76.30 1.58 
10 90 30 ee 1.60 
10 90 30 79.21 1.61 
10 90 30 78.46 1.58 
10 90 30 76.68 1.55 
10 90 30 74.75 1.63 
10 | 90 30 74.71 1.64 
30 90 30 134.15 1.52 
30 90 30 136.34 1.51 
30 90 30 135.84 1.50 
30 90 30 137.40 1.53 
30 90 30 137.34 1.50 
30 90 30 135.24 1.48 
30 90 30 133.03 156 
30 90 30 | 133.62 1.52 
30 90 30 130.54 1.54 
30 90 30 132.37 | 1.51 
50 90 30 193.03 1.39 
50 90 30 194.05 1.48 
50 90 30 187.17 1.50 
50 90 30 174.42 1.49 
50 90 30 187.96 1.39 
50 90 30 188.53 1.48 
50 90 30 180.46 1.51 
50 90 30 184.40 1.50 
50 90 30 182.35 1.47 
50 90 30 179.00 1.48 
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EVAPORATION TIME AND SHAPE FACTOR OF THE DROPLETS WHICH HAVE 
BEEN SELECTED AS REPRESENTATIVE CASES © 


Initial Ternperature Evaporation Time Shape Factor B 
of the solid surface (°C) (s) (adirn.) 
50 microliters 
180 14.4 3.18 
164 Boe i 2.76 
155 453 2.39 
150 42.5 2.34 
145 52.8 - 2.46 
130 86.2 1.82 
110 110.8 171 
90 194.1 1.48 
30 microliters 
180 Laez 3.34 
164 31.0 2.12 
55 32.9 252, 
150 35.4 2.38 
145 41.4 2.30 
130 58.9 1.93 
110 85.7 Lee 
90 137.4 iS) 
10 microliters 
180 8.4 3.39 
164 15 2.94 
155 16.5 3.08 
150 20.2 2.60 
145 723) Peay 
130 22 2.05 
110 47.4 1.80 
90 79.2 . 1.61 
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INITIAL TEMPERATURE OF THE SOLID SURFACE, OF THE CHILL PLATE AND 
OF THE HEATERS FOR EACH TEST SERIES 


Initial Temperature Chill Plate Temperature 
of the solid surface (°C) Temperature (°C) of the Heaters (°C) 
90 30 522 
110 33 635 
130 34 655 
145 35 675 
150 35 702 
155 36 715 
164 | eH 730 
180 59 750 
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PERCENTAGE OF HEAT FLUX DIRECTLY ABSORBED BY THE DROPLET 


10 microliters (vaporization heat = 22 joules) 


HEATERS 
TEMP. 
(°C) 


750 
730 
il 
702 
675 
655 
635 
Sys) 


~ SURF. 
TEMP. 


(°C) 


180 
164 
155 
150 
145 
130 
110 
90 


INCIDENT 
FLUX 
(KW/M") 


att 
10.2 
9.6 


aol 


8.2 
TE 
6.9 
4.1 


30 microliters (vaporization heat = 67 joules) 


HEATERS 


SURF. 


INCIDENT 
FLUX 


(KW/M?) 


ie 
10.2 
9.6 
val 
8.2 
ie 
6.9 
4.1 


BETA 


3.39 
2.94 
3.08 
2.60 
Zh 
2.05 
1.80 
1.61 


BETA 


3.34 
DHT? 
2.52 


.2.38 


2.30 
19S 
1.71 
bin 


ABSORBED %OF 
HEAT 


HEAT 
(J) 


oy! 
6.5 
72 
6.0 
5.9 
4.8 
5.0 
3.8 


ABSORBED 
HEAT 


J) 


13.8 
23.4 
20.0 
18.2 
17.6 
16.1 
16.8 
Zz 


50 microliters (vaporization heat = 112 joules) 


HEATERS SURF. | INCIDENT EVAP. 


TEMP. TEMP. FLUX TIME 
(°C) (°C) (KW/M*) (s) 

750 180 11.1 14.4 
730 164 10.2 33.1 
715 155 96 453 
702 ater bof 9.1 42.5 
675 145 8.2 52.8 
655 130 - TS 86.2 
635 110: 6.9 - 110.8 
525 90 4.1 194.1 


BETA 


3.18 
2.76 
aso 
2.34 
2.46 
1.82 
Wert 
1.48 


ABSORBED 
HEAT 
(J) 


nice 


28.0 
24.8 
214 
2571 
21.0 
Pikes 
16.3 


Jo 


OF 


HEAT 
ABS. 
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ABSORPTION COEFFICIENT VALUES FOR WATER IN THE RANGE 1-10 MICRON 
(Siegel et al. 1981) 


Wavelength (micron) Absorption Coefficient (cm) 
1.0 3.63 x 107 
12 1.04 
1.4 1.24 x 10! 
1.6 6.72 
1.8 8.03 
2.0 . 6.91 x 10! 
we) 1.65 x 10! 
2.4 5.01 x 10! 
2.6 1.53 x 10? 
2.8 | | 5.16 x 10° 
3.0 Tidexel0: 
3.2 i 3.63 x 10° 
S34 ae fpAls 2itr 
3.6 1.80 x 10? 
3.8 1.12) xe 
4.0 1.45 x 10? 
4.2 | 2.06 x 10? 
4.4 2.94 x 10? 
4.6 4.02 x 10° 
4.8 7 3.93 x 10? 
5.0 3.12 x 10? 
55 2.65 x 10° 
6.0 224 x 10° 
6.5 758 x 10? 
7.0  wbe/4 102 
US 5.46 x 10° 
8.0 5.39 x 10? 
8.5 5.43 x 10? 
9.0 5.57 x 10? 
9.5 5.87 x 10? 
10.0 6.38 x 10? 
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RADIATION ABSORPTION IN WATER 


The water droplets absorb most of the radiation present on the surface of the water 
droplets. The following is a chart showing the percentage of radiation which is absorbed 
by water for various thicknesses and various spectrums produced by the heaters at different 
temperatures.. These numbers were calculated by a simple computer program. 


THICKNESS HEATERS © RADIATION 
(mm) | TEMP. (°C) ABSORBED (%) 
1.0 750 97.4 
1.0 w20 97.7 
he) 700 98.0 
1.0 675 98.2 
1.0 650 98.4 
1.0 625 98.7 
1.0 600 98.8 
0.5 750 94.8 
oS q25 95.3 
0.5 700 95.8 
0.5 675 96.3 
0.5 650 96.7 
0.5 625 97.1 
0.5 600 97.4 
0.1 750 82.0 
0.1 10) 83.0 
0.1 700 84.0 
0.1 675 85.0 
0.1 650 86.0 
0.1 625 86.9 
0.1 600 87.8 
0.01 750 413 
0.01 tps 41.7 
0.01 700 42.1 
0.01 675 42.5 
0.01 650 42.9 
0.01 625 43.2 
0.01 600 43.5 
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PART C 


MULTI-DROPLET MODEL 


- 1 


Cl. INTRODUCTION 


The spray cooling associated with the solid fuel fire extinguishment process is numerically 
simulated in this portion of the report. On the basis of the results provided by the single- 
droplet code, a data base is built to extend the analysis to a multi-droplet scenario. A 
closed form solution is developed to approximate the surface temperature at times and 
locations remote with respect to the point considered. 

The surface temperature calculation in the multi-droplet model is performed by 
superimposing the appropriate elements selected from the data base. A detailed 
temperature distribution on the solid surface under the spray cooling is obtained. 

A spatial average temperature of the solid surface is calculated and its significance is 
discussed in relation with the influence of dominating parameters such as droplets mass flow 


rate, impingement area and initial solid surface temperature. 


C2. MULTI-DROPLET MODELING 


In the study of spray cooling, a major objective is the determination of the average 
temperature of the solid surface as a function of time. Eventually, one would alo need to 
know when a steady state temperature is reached for a given spray mass flow rate. Other 
information of interest concern the solid surface temperature distribution, the initial cooling 
effect of the impinging spray, etc. 

In order to model a multi-droplet evaporative cooling transient, two basic steps have to be 


completed: the droplet distribution simulation and the temperature distribution calculation. 


C2.1. Droplet distribution simulation 

As well as for single droplet cooling, also for multi-droplet situations the solid surface starts 
from an initial condition of uniform temperature. During the cooling process, Te droplets 
are randomly deposited (in terms of landing time and position) on the surface. In order to 
obtain a suitable droplet distribution pattern for later use in multi-droplet temperature 
calculations, a numerical model is created which simulates a random droplet distributer. 
The program is able to provide the information such as the deposition time and the location 


on the surface for all droplets. 


The distributer moves with initial constant speed within the boundaries of a circle of fixed 


radius, hits the walls and bounces back on a new trajectory like in a "Pin-ball" game (see 


Cea 


Fig. C1). During this movement, it releases the droplets with constant frequency. The 
numerical simulation of this behavior can be described as follows. The droplet dispenser 
is represented by a moving point, whose initial position and initial velocity are assigned. It 
starts to move freely inside the circle and delivers a droplet at each given time interval. 
This behavior 4 simulated by recording the point location at the moment corresponding to 
the fall of each droplet. During each time period, the trace of the point is followed by 
tracking the position and velocity of the point in small time intervals. Periodically the point 
approaches the boundary. Assuming that all contacts are elastic, the rebounding angle and 
velocity can be computed by applying the principle that the nae angle is equal to the 
rebounding angle and that the velocity remains unchanged. It has to be pointed out that 
when the point passes through the center of the Bea: or very close to that region, it will 
keep rebounding within that area, since the impact and rebounding angles become very 
small. In order to avoid such repeating moves, one can force the point velocity vector to 
change its direction when it happens to pass through the center point of the circle or the 
nearby region. Further improvement can be obtained by assuming that the circle boundary 
wall rotates against its center point with a time dependent frequency, for example a 
sinusoidal function, so that the rebounding angle can be different from the impact angle. 
An additional feature is included in the model in order to avoid overlapping among the 
droplets which have been deposited on the same or very nearby location. This is due to the 
purpose of applying the results obtained from the single droplet model, as it will be 
discussed later. The program simply forces the distance between those droplets which are 


deposited within a given period to be no less than a twice the size of a single droplet 


Ga 


diameter. 
A typical plot showing the random distribution of droplets on the surface which is obtained 


with this numerical simulation is displayed in Fig. C2. 


C2.2. entrees tire distribution calculation 

The solid surface which is interested by the cooling effect of the spray evaporation is the 

same that is described for the single droplet cases, i.e., its initial temperature is constant and 

uniform, and its depth and width are big enough to be considered as infinite. The other 

major assumptions used in both single- and multi-droplet modeling are: 

1) The wetted area between each liquid droplet and the solid surface is a circle and 
remains constant during the evaporation process. 

2) At any moment, each droplet has a spherical segment shape. 

3) Nucleate boiling is suppressed in the liquid-solid interface. 

4) Conduction is the only effective heat transfer mode considered inside the liquid 
droplets. 

Since the paneiction is the dominant heat transfer mode, the governing equation for this 


multi-droplet cooling can be written, for both the droplets and the solid, as: 


OT 
Lee oe er 
oe aV°T (C1) 


Similarly to the single droplet model, the initial conditions can be expressed as: 
The T ater = 0 (C2) 
For each droplet, the corresponding initial surface temperature is different because of the 
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different time and location of impingement. This condition can be written as: 


T = T; at “loca lone pee Cepia a eaten (C35 


+30 
This condition is used during the entire multi-droplet evaporation modeling, as discussed in 
detail in later sections. The assumption that the surface temperature T,,, in each droplet 
region at the initial moment is uniform and equal to the temperature at the center point of 
the wetted area is implicitly used. This assumption is reasonable if the wetted area under 
each droplet is small, as well as the temperature gradient in this region. From the surface 
temperature profiles in the single droplet model, it can be seen that the large temperature 


gradients on the solid surface occur near the droplet edge. 


In analogy with the single droplet case, the boundary conditions are assigned as: 


e =0 when Ir-© | (C4) 

= = 0 when z-7-© (C5) 

Kal palh =k, er at Z=0 droplets locations (C6) 
-k, 2 = Tiel) at z=0 elsewhere (C7) 


In general, the temperature difference used here is defined as: 


where T,, is the intial solid surface temperature. By using Eq. (C8), the governing equation 
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Ue LT y0 (C8) 


and the initial and boundary conditions become: 


OU ss 2 
55 a V°U (C9) 
T= 0 aly enb =i) (C10) 
OU _ is 
aCe 0 when (G11) 
OU _ cag 5 
3s 0 when 2Z (C25 
OU oU A 
ky (=) -=X; (==), at 2=0, droplets locations (C13) 
Oz dz 
-c,2U = nu at 2=0 elsewh 
gage at z=0: jetsewnere (C14) 


The governing equation and the associated conditions can be solved for the droplets and the 
solid surface temperature by using the Boundary Element Method (BEM) and the Control 
Volume Method (CVM) as discussed in the single droplet modeling presented in Part A and 
more extensively in a previous report (Tartarini et al., 1990). Unlike the single droplet 
calculation, this three dimensional conduction problem can not be reduced to a two 
dimensional case by applying the symmetry for both the droplet and the solid surface. For 


thousands droplets on the relatively large solid surface, the calculation procedure would be 
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huge in terms of computer memory and CPU. 


The governing conduction equation (C1) applied to a multi-droplet scenario constitutes a linear 


system; therefore, superposition can be used for this multi-droplet modeling. 


Assume JU; is the temperature difference for a single droplet , deposited on a solid surface 


at uniform initial surface temperature T.,. at the time ¢;: 


The governing equation for the solid in single droplet model is written as: 


OU, 
=a = «V*U; 
and the associated boundary conditions are: 

OU; 

AG = 0 when r- 

OU; 

= = 0 when z-© 

OU; OU; 

kK aot muki AE atiyjez = 
OU, 
a = AU; elsewhere 


The initial condition is: 


(C15) 


(C16) 


(C17) 
(C18) 
(C19) 


(C20) 
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U; = 0 at t=t,;, i=1,ND voer) 
Now one can assume: 
ND 
v=) y, | (C22) 
1=1 


Thus, the governing equation becomes: 


OU. 2 
aa Civ - Lh (C23) 
and the boundary conditions are: 
oun sit 
epee 0 when (C24) 
Omt> WO) 
Seen 0 when 2Z CCZ5) 
k (52) =k (32 at droplets locations (C26) 
pyoc me \ Oz), 
alts hu elsewhere (C27) 
oz 
The initial condition is: 
U=0 at. t=0 (C28) 


These equations coincide with those which have been written for the multi-droplet model. 


Therefore U, as it is defined in Eq. (C29), is indeed the solution for the multi-droplet 


cooling on the solid surface: 


U =" Toit 5 = eae) ) el Tia) (C29) 


or, in more general form: 


ND 
dese Gila Ya i acter i (C30) 


i=1 


It can also be shown that the heat flux for the multi-droplet model can be calculated by 


superposition due to the linear relation between the temperature and the heat flux. 


There are different approaches for solving the multi-droplet problem by applying the 
superposition technique. However, since the studies presented in this section of this report 


are still in progress, only a summary of the procedures utilized until now is contained here. 


The single-droplet code has been used to calculate transient temperatures and heat fluxes. 
A data base has been built with the results from the single-droplet code. Then, in order to 
avoid the huge computational expenses connected with the run of the single-droplet code 
for thousands of droplets, proper values for the temperatures at various time step and 


location have been selected from the data base. 


In order to complete the data base, one has to calculate the temperature difference 7-7, 


for those locations and times which are not contained in the part of the data base originated 
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by the single-droplet code. This is obtained by means of closed form solutions which had 
to be created for this specific purpose. For the present purposes, it is only necessary to have 
the final expressions of those solutions, which are based on the following expression for the 


temperature difference for a semi-infinite solid surface: 


V Py; A 
—— 
(es 2 
GT pork Tas ~__Ps Sps_oyp{-_2 (C31) 
a Sr oe 4a,t 


The first solution is the recovery function, which yields the temperature difference 7-7, for 
the far-field region after the complete evaporation of the first droplets. This function is 


obtained by substituting a modified time t-1 into Eq. (C31). The resulting expression is: 


p> A 
es Ge at Vas, tee e i r? (oeoe) 
ee ae ice ee mA aD 


The second solution is the evaporation function, which yields the temperature difference in 
the far-field during the evaporation process. It is obtained as a modification of the recovery 


solution. By using a small value of time ét, Eq. (C31) can be written as follows: 


P, A 
Lrg C r? T t (C33) 
= = pA phan vs ter DS ee aspen tee my PS pei 
oo a | re ez}o9($(2 4) 


where the factor 


(3-4) 


is 0 at t=O and 1 at t=7. 


Figure C3 shows a comparison between the calculated temperature and the closed form 


solutions. 


The full data base consists of results from the single-droplet model and calculations 
performed with the closed form solutions mentioned above. At this point, the data base can 
provide all the required information for applying the superposition criterion to the multi- 
droplet modeling. The major advantage related to the use of this data base is the saving 
of computational time. Since all the single-droplet calculations are performed at a 
preliminary stage, the execution of the multi-droplet program is very rapid. On the other 
hand, this results in the need of a very large computer storage where the data base has to 


be contained. 
The multi-droplet program, which has been created on the basis of the concepts above 


summarized, is listed in the section C6. The latest version of the single-droplet code is also 


included. 
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C3. RESULTS 


The main objective of the multi-droplet model is to calculate the temperature distribution 
of the solid surface as a function of time during a spray evaporative cooling. As described 
in the previous section, the temperature at any point of the solid surface can be obtained 


by superimposing the results from many runs of the single-droplet code. 


Given an initial solid surface temperature and a pre-determined droplet distribution pattern, 
the detailed temperature distribution at any selected time can be obtained. Figures C4, C5 
and C6 show some typical three dimensional temperature plots for a portion of the solid 
surface at different times. 

During the first seconds of the transient, the region influenced by the cooling effect of a 
single droplet is easily individuated and, at a sufficient distance from the droplet, the solid 
temperature remains constant and uniform. Consequently, the temperature drop is very 
sharp in the areas where a droplet is evaporating. After few seconds, more droplets are 
deposited on the solid surface, and the temperature at any point becomes a combination of 
the cooling effect generated by all the droplets together. At this point, the single droplets 
influence cannot be individuated as easily as before. After a longer time, the plot shows an 
overall temperature drop and the solid surface temperature is everywhere below its initial 
value. 

Obviously, the lowest solid temperature at any time is always higher than the initial 


temperature of the liquid droplets (20°C). 


G13 


The most useful parameter to describe the temperature distribution for multi-droplet cooling 
is the spatial average temperature of the considered portion of the solid surface. It is 


defined as follows: 


Npi Np; 


[ {taxay DD Tig Axi Ay; 

= — i2=1 j=1 

Teyg = = (C34) 
ff dxdy 


yy e. Ax,Ay, 


pital SEAL 


Further simplification can be obtained if the same distance is imposed between the grid 
points in both directions of the mesh. Then the spatial average temperature is interpreted 
as the temperature summation for all the considered points on the surface divided by the 


total number of guid points: 


Np 

de ok (C35) 
Tague fe 
avg N 


Figure C7 shows an average temperature curve versus time. In this plot, the average 
temperature is calculated at every second. Fluctuations of the average temperature are 


observed. 


Since the droplets distribution is irregular in time and space, the cooling effect around a 
single point cannct be monotonic. This explains the oscillatory behavior of the average 


temperature curve, that could be smooth only if the droplets distribution were uniform both 
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in time and in space. 


The governing parameters for the multi-droplet cooling include the initial solid surface 
temperature, the physical properties of the solid surface and the mass flow rate of the 
droplets. In this study, only Macor is considered, but other materials characteristics are 
already implemented in the code for future tests. The initial solid surface temperature is 
always less than 160°C to prevent boiling, as described in the previous sections, i.e. Part A 
and B. The mass flow rate is added to the cooling eno which refer to the properties 
of the solid. Without losing generality, one can assume that all the droplets have the same 
initial volume Vp. Under this assumption, the definition of the mass flow rate can be 


written as: 


(C36) 


where f is the droplets depositing frequency, and A is the area over which the averaging 
process is performed. In this definition, the liquid density and the droplet initial volume are 
assumed to be constant values, so that the effective variables for the mass flow rate are the 


frequency and the area. 


Figures C8 to C10 show some different average temperature curves for a common initial 
solid temperature of 160°C, a constant total area, A1=0.0314 m’, and different frequencies 
of 1, 2, and 5 droplets per second respectively. Figures C11 to C13 show other average 


temperature curves for A2 = 2 Al, and a frequency of 2, 4, and 10 droplets per second 
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respectively. Figures C14 to C16 display the curves for A3 = 4 Al and frequencies of 4, 8 
and 20 droplets/s. The higher is the frequency, the lower is the average temperature. 


If one defines the droplet distribution density is: 


52) 


then the mass flow rate can be expressed as: 
h = pV,d (C38) 


Figure C17 illustrates some different average temperature curves versus the number of 
droplets. For this plot, the droplet distribution density, as well as the mass flow rate, is 
constant. The mass flow rate appears to be a dominant factor, since all the average 


temperature plots show the same trend for the same value of the mass flow rate. 


All the average temperature plots indicate that the average temperature tends to reach 
constant values, 1.e. a steady state average temperature. The existence of the steady state 
average temperature can be demonstrated by applying an energy balance to a control 
volume of the solid surface of concern. From the first law of thermodynamics, the energy 


balance is expressed as: 


90 ca cuaeng 
ot Qin Qout i 


where Q is the total heat stored in the control volume, which can be written as: 


Q =f Peet av (C40) 
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By selecting very small and constant thickness dz for the control volume, being dV = dz dA, 


the expression for the total heat becomes: 


O = pcpdz [ TaA (C41) 


Recalling the definition of the average temperature, the total heat Q can be further 


expressed as: 


Ore pce dzan mia. (C42) 


so that the energy balance equation (C46) becomes: 


Tai 
(p €, dz Ay See an O TF - Dore (C43) 


When the outgoing heat flux due to the evaporation is equal to the incoming heat flux 
supplied from the rest of the semi-infinite solid, the right hand side of the above equation 


becomes zero, and one has: 


SSZVENO (C44) 


Then a steady state average temperature is reached. 


Figures C18 to C20 show some average temperature curves for the same mass flow rate at 
different initial solid temperatures. They all show a quick temperature decrease in initial 
portion of the transient. This can be explained if one remembers that the evaporation time 
for Macor under these initial conditions is in the range between 25 and 100 seconds; 


therefore, during this initial period the solid surface is consistently cooled by the evaporating 


(Oz558 | 


droplets, and the average temperature continuously decreases. After this, some of the 
droplets are completely evaporated and the recovery process starts in those area previously 
covered by the water, thus the temperature increases again locally. Since there are always 
some evaporating droplets and recovery areas in the same region, the combined process 
results in a relatively smooth and slow variation in the average temperature. The steady 


State average temperature is approached as time tends to infinity. 
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C4. CONCLUSIONS 


A multi-droplet model has been developed to obtain the transient temperature distribution 
on a hot solid surface subjected to evaporative cooling. The numerical simulation can be 
applied to a broad range of solid materials, from low thermal conductivity materials like 


Macor to high thermal conductivity materials like aluminum. 


The single-droplet model, which has been previously validated, provides numerical 
predictions of the evaporation time. It also provides temperature distribution, as well as 
heat flux, inside the droplets and along the solid surface. These results are used to generate 


the data based used by the multi-droplet program. 


By applying the characteristic of linearity of the conduction equation and assuming 
homogeneous boundary conditions along the exposed solid surface, the superposition 
technique is used for the multi-droplet simulation. Given a pre-determined: droplet 
distribution pattern on the solid surface, the temperature distribution is obtained by 
superposition of results from the single-droplet model with different initial and boundary 
conditions. The summation of these contributions provides the detailed temperature 


distributions as a function of time during the multi-droplet evaporative cooling. 


A spatial average temperature is used to summarize the cooling effect for the solid surface. 


The influence of the spray mass flow rate and the initial solid temperature is studied and 
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it is observed that the mass flow rate is a dominating factor in the multi-droplet evaporation 


process. A steady state value of the spatial average temperature is also calculated. 
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C5. FIGURES 
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Fig. C2 Typical distribution of droplets on the surface 
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Fig. C4 Temperature distribution on the solid surface; t = 10s, f = 2 droplets/s, 
Vy =5 ul 
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Fig. CS Temperature distribution on the solid surface; t = 10s, f = 10 droplets/s, 
Vo = 5 wl 
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C6. APPENDIX - MULTILDROPLET CODE 


In the following pages, the current version of the multi-droplet code is presented. It 
includes: a) the main program for multi-droplet calculations; b) the program calculating the 
recovery process which follows the evaporation of a droplet and c) the single-droplet 


program which is used to build the multi-droplet data base. 
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- -MULTI-DROPLET MODEL (MACOR) 


DIMENSION U(41,41) ,TR(3000,5) 
DIMENSION UDATA(27,1,1,15,973) , TEND(27) , TEVAP (27) 


DIMENSION PT(50) 


. -NIX=41,NIY=41,NDROP=30, MNTSO=20,MNVO=1 , MNBETA=1 , MNRL=30 , MNTL=60 


...DIMENSION U(NIX,NIY) , TR(NDROP,5) , TEND(MNTSO) , TEVAP (MNTSO) 
. ..DIMENSION UDATA(MNTSO,MNVO,MNBETA, MNRL,MNTL) 

COMMON /DDATA/DTSO,DTL, DRL, TO, DTIME 

COMMON /DROP/TSO,VO,BETA,TL 

COMMON /NDATA/MNTSO,MNVO,MNBETA, MNRL, MNTL 


COMMON /NUMR/NRL,NOD,NDROP 


COMMON /GEOM/BLX, BLY, OBX,OBY,DBLX, DBLY 


COMMON /XY/NIX,NIY,NOX, NOY 


COMMON /FR/XMO,YMO,TDROP, FLOWR 

OPEN (UNIT=66, FILE='MD1. DAT! ,STATUS='NEW') 
C..GIVE INITIAL SOLID SURFACE TEMP. TO 

PRINT’ *, "INITIAL SOLID TEMP (< 160 C) TO=?' 


READ(5,*) TO 

.PLOT INFORMATION 
PRINT *,'HOW MANY PLOTS?! 
READ (5,*)NNPT 


PRING: *40° PLOT AT (SECOND) "2s 


READ(5,*) (PT(I) ,I=1,NNPT) 


.NUMBER OF DROPS PER SECOND (STANDARD 10 D/SEC) 


NDS=10 


PRINT *,'NUMBER OF DROPS PER SECOND NDS='! 


READ(5,*) NDS 
FNDS=10./(1.*NDS) 
. GEOMETRY 
..LOWER LEFT POINT OF THE BLOCK 
OBX=121.4214 
OBY=121.4214 
OBX=80. 
OBY=80. 
..SIZE OF THE BLOCK 
BLX=40. 
BLY=40. 
..SIZE OF THE CELL 
DBLX=1. 
DBLY=1. 


NDROP=3000 
NOD=0 


MNTSO=27 

MNVO=1 

MNBETA=1 
- VO=5 

MNRI=15 
- VO=5 

MNTL=973 
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KPT=1 
Cu 

CALL MDXY ; 
c 2m 

CALL MDTR(TR,FNDS) 
Cae 

CALL MDDATA(UDATA, TEND, TEVAP) 
Gai. 2 


DO 15 ID=1,NDROP 
TR(ID,2)=TO 
15 CONTINUE 


Cae. 
TIME=0. 
DTIME=20. 
TSTOP=(1.*NDROP) /(1.*NDS) 
Ca 


C...MASS FLOW RATE CALCULATION (KG/M‘2 S) 
PI=4.*ATAN(1.) 
RHOW=996. 
AREA=PI* (XMO*1.E-3) **2 
VOLUME=NDROP*V0*1.E-9 
TMASS=VOLUME*RHOW 
FLOWR=TMASS/TDROP/AREA 


PRINT *,'NUMBER OF DROPLETS=',NDROP,' VO0=',VO,'(MM‘3) '! 
$ ' TOTAL TIME=',TDROP,'(SECOND) ' 
PRINT *,'THE MASS FLOW RATE IS',FLOWR,'(KG/M‘2 S)'! 
eli 
1 CONTINUE 
cue 


CALL MDIN(TIME,TR,U,UDATA, UMD, TEND, TEVAP) 
C...SPACIAL AVERAGED TEMP. IN BLOCK 
TAVE=0. 
DO 50 J=1,NIY-1 
DO 50 I=1,NIX-1 
TAVE=TAVE+ . 25* (U(I,J)+U(I+1,J)+U(I,Jd+1)+U(I+1,J+1) )+T0 
50 CONTINUE 
TAVE=TAVE/ ( (NIX-1) * (NIY-1) ) 
PRINT *, 'TAVE=',TAVE 
WRITE (66, *) TIME, TAVE 
Cue PLor 
NPT=0 
DO 70 IPT=1,NNPT 
IF (PT(IPT) .EQ. TIME) NPT=1 
70 CONTINUE 
IF (NPT.EQ.1) THEN 
PRINT *,'A PLOT AT TIME=',TIME 
CALL MDPLOT (KPT, TIME, TAVE,U) 
ENDIF 


LE (TIME EO. STOPyGCOsLO m2 


TIME=TIME+DTIME 


GOsTO «1 
2 STOP 
END 
C....SUB. FOR DEFINING THE GRID POINTS IN THE TEMP. CALCULATING 
BLOCK FOR 
C....-PROGRAM MD.FOR 


SUBROUTINE MDXY 
COMMON /DDATA/DTSO,DTL,DRL,TO,DTIME 
COMMON /DROP/TSO,VO,BETA, TL 
COMMON /NDATA/MNTSO,MNVO,MNBETA,MNRL, MNTL 
COMMON /NUMR/NRL,NOD,NDROP 
COMMON /GEOM/BLX, BLY, OBX, OBY, DBLX, DBLY 
COMMON /XY/NIX,NIY,NOX,NOY 

C...INDEX NUMBERS FOR CORRESPONDING X AND Y DIRECTION 
NIX=INT (BLX/DBLX+1+.000001) 
NIY=INT (BLY/DBLY+1+.000001) 
NOX=INT (OBX/DBLX+1+.000001) 
NOY=INT (OBY/DBLY+1+.000001) 


RETURN 

END 
C...-SUB. FOR ASSIGNING THE DROPLETS DISTRIBUTION IN TERMS OF THE 
TIME THEY 


C....ARE DEPOSITED ON THE SOLID AND THE LOCATION 
C..TR(ID,1):TIME TO START 
C..TR(ID,2):INITIAL SURFACE TEMP 
C..TR(ID,3):X 
C..TR(ID,4):Y 
C..TR(ID,5):INDICATOR OF INSIDE OR OUTSIDE THE CALCULATION BLOCK 
SUBROUTINE MDTR(TR, FNDS) 
DIMENSION TR(NDROP, 5) 
COMMON /DDATA/DTSO,DTL, DRL,TO,DTIME 
COMMON /DROP/TSO,VO0O,BETA, TL 
COMMON /NDATA/MNTSO,MNVO,MNBETA,MNRL, MNTL 
COMMON /NUMR/NRL,NOD,NDROP 
COMMON /GEOM/BLX, BLY,OBX,OBY, DBLX, DBLY 
COMMON /XY/NIX,NIY,NOX,NOY 
COMMON /FR/XMO,YMO,TDROP, FLOWR 
C...MDTR10.DAT: R=100 MM 
C...MDTR2A.DAT: R=141.4124 MM 
C...MDTR200.DAT: R=200 MM 
OPEN (UNIT=15, FILE='MDTR2A.DAT' , STATUS='OLD') 
C...DISTANCE FROM CENTER OF THE PLATE TO THE EDGES IN X AND Y DIR. 
READ (15, *) XMO, YMO 
C...XMO,YMO (MM) 
READ (15, *)NNDROP 
DO 10 I=1,NDROP 
READ(15,*) ID 
B... + LD=T 
READ(15,*) TID 
READ(15,*)XI,YI 
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TR(I,1)=TID*FNDS 
C...TR(I,2):TEMPIDO IS NEEDED TO BE DETERMINED. 
C....MM 

TR(I,3) =XI+XMO 

TR(I,4)=YI+YMO 

TR(I,5)=0. 
C...IO=1:INSIDE; IO=0:OUTSIDE 

IF((TR(I,3).GE.OBX.AND.TR(I,3) .LE. (OBX+BLX) ) .-AND. 

$ (TR(I,4).GE.OBY.AND.TR(I,4).LE. (OBY+BLY) ) )TR(I,5)=1. 
10 CONTINUE 

TDROP=TR(NDROP, 1) -TR(1,1) 

CLOSE (15) 

RETURN 

END 
C...SUB. FOR CREATING A DATABASE FOR TEMP. DISTRIBUTION 
C...AT DIFFERENT INITIAL SOLID SURFACE TEMP. (TSO), 


Gace INITIAL DROPLET VOLUME (VO) 
Oe BETA (BETA) 
Gut LOCAL TIME (TL) 
care LOCAL RADIUS (RL) 
SUBROUTINE MDDATA(UDATA, TEND, TEVAP) 
G PARAMETER (MNTSO=6,MNVO=1,MNBETA=1,MNRI=4 /MNTL=4) 


DIMENSION UDATA(MNTSO,MNVO,MNBETA , MNRL, MNTL) 
DIMENSION TEND(MNTSO) , TEVAP(MNTSO) 
COMMON /DDATA/DTSO,DTL, DRL, TO, DTIME 
COMMON /DROP/TSO,VO,BETA,TL 
COMMON /NDATA/MNTSO,MNVO,MNBETA,MNRL, MNTL 
COMMON /NUMR/NRL,NOD,NDROP 
COMMON /GEOM/BLX, BLY, OBX,OBY,DBLX, DBLY 
COMMON /XY/NIX,NIY,NOX,NOY 

C...MDDATA2.DAT: VO=5, BETA=1.4,TSO=160 - 30 
OPEN (UNIT=25, FILE='MDDATA2.DAT!', STATUS='OLD') 


Chee UDATA: TEMP. (C) 
C......ITSO=1: TSO=160 (C) 
Cheer 2: TSO=155 

eis Sais 3: TSO=150 

Cc, Be. NVO=1: VO=10 (MM‘3) 
Cun eee NBETA=1: BETA=1.4 
C......RL=NRL*DRL (MM) 
C......TL=NTL*DTL (SECOND) 
CM 4 25. ...'M': MAXIMUM 
oe 


C...FOR IVO=1 
READ(25,*) VO 
READ (25,*) TTO 
READ (25,*) DTSO 
READ (25, *) DTL 
READ(25,*) DRL 


READ (25, *)NTSO 
IVO=1 
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corn 


IBETA=1 
C...FOR IBETA=1 

BETA=1. 4 
Cc. 

READ (25,*) NRL 

DO 10 ITSO=1,NTSO 

READ (25, *) IITSO 

READ (25, *) TSOI 

READ (25, *)NTEVAP, NTEND 

TEVAP (ITSO) =NTEVAP 

TEND (ITSO) =NTEND 
Cee: 


DO 20 ITL=1,NTEND 
READ (25, *) (UDATA(ITSO,IVO,IBETA,IRL,ITL) , IRL=1,NRL) 
DO 21 IIRL=1,NRL 
IF (UDATA(ITSO,IVO,IBETA, IIRL,ITL) .GT.0.) THEN 
UDATA(ITSO,1IVO,IBETA, IIRL, ITL) =0. 
ENDIF 
21 CONTINUE 
20 CONTINUE 
ahr 
10 CONTINUE 
CLOSE(25) 
RETURN 
END 
C...SUB. FOR CALCULATING THE TEMP. ON ALL GRID POINTS IN THE 
CALCULATION 
C. seBLOCK: 
SUBROUTINE MDIN(TIME,TR,U,UDATA, UMD, TEND, TEVAP) 
DIMENSION U(NIX,NIY) , TR(NDROP, 5) 
DIMENSION UDATA(MNTSO,MNVO,MNBETA,MNRL, MNTL) 
DIMENSION TEND (MNTSO) , TEVAP(MNTSO) 
COMMON /DDATA/DTSO,DTL, DRL, TO, DTIME 
COMMON /DROP/TSO,VO, BETA, TL 
COMMON /NDATA/MNTSO,MNVO,MNBETA,MNRL,MNTL 
COMMON /NUMR/NRL,NOD,NDROP | 
COMMON /GEOM/BLX, BLY, OBX, OBY, DBLX, DBLY 
COMMON /XY/NIX,NIY,NOX,NOY 
EXTERNAL UFEV 
EXTERNAL UFRE 
DO 11 I=1,NIX 
DO 11 J=1,NIY 


U(I,J)=0. 
11 CONTINUE 
TSO=TO 
Cc PRINT *,' NOD=',NOD,' TR(ID,2):! 
Cc PRINT *,(TR(ID,2) ,ID=1,NOD) 


DO 10 ID=1,NOD 
C...DETERMINE TSO,VO,BETA,TL 
TSO=TR(ID,2) 
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IF(TSO.LT.30.)PRINT *,'TSO=',TSO, 
$ 'MORE DATA NEEDED FOR TSO BELOW 30 C'! 
KTSO=0 
DO 12 K=1,MNTSO+1 
IF (TSO.LE.TO-(K-1) *DTSO.AND.TSO.GT.TO-K*DTSO) KTSO=K 
12 CONTINUE 
C...TL=(ITL-1) *DTL 
TL=TIME-TR(ID,1) 
C..FOR EVERY POINT IN THE BLOCK, DETERMINE RL 
DO 20 I=1,NIX 
DO 20 J=1,NIY 
RL=SORT( (TR(ID, 3) —(NOX+I-1) *DBLX) **2 
+(TR(ID,4) -(NOY+J-1) *DBLY) **2) 
C...RL=(IRL-1) *DRL 
e IF (RL. LE. (MNRL-1) *DRL.AND.TL. LE. MNTL*DTL) THEN 
IF (RL. LE. (MNRL-1) *DRL.AND. 
S TL. LE. TEND(KTSO) .AND. TL. LE. TEND(KTSO+1) ) THEN 
CALL MDU(UDATA,TSO,TL,RL,UMD) 
UIJ=UMD 
ELSE 
C..CALL FUNCTION 
IF (TL. LT. TEVAP(KTSO) ) THEN 
UIJ=UFEV(TSO,V0,BETA, TL, RL, KTSO, TEVAP, TEND) 
ELSE 
UIJ=UFRE(TSO,V0,BETA,TL,RL,KTSO, TEVAP, TEND) 
ENDIF 
ENDIF 
C.. SUMMATION 
U(1I,J)=U(I,J)+UIT 
20 CONTINUE 
10 CONTINUE 
C..DETERMINE TR(ID,2) AND NEW NOD (NUMBER OF ON-RUN DROPLETS) 
NOD1=NOD 
TSO=TO 
DO 100 ID=NOD1,NDROP 
UOUT=0. 
C..WHEN THE DROPLET IS DEPOSITED 
IF(TR(ID,1).GE.TIME.AND.TR(ID,1) .LT.TIME+DTIME) THEN 
C..PTS INSIDE OF BLOCK 
IF(TR(ID,5) .EQ.1) THEN 
I=INT(TR(ID,3) /DBLX) -NOX+1 
J=INT(TR(ID,4) /DBLY) -NOY+1 
C..INTERPOLATION TO GET TR(ID,2) 


C....FROMsU Gly J) PAU(LT,J+1) ,U(I+1,3+1) AND U(I+1,d) 
UX1=(U(I+1,J7)-U(I,J) ) *(TR(ID,3) -(NOX+I-1) *DBLX) /DBLX 
$ +U(I,d) 


UX2=(U(I+1,J+1) -U(I,J+1) )*(TR(ID,3) -(NOX+I-1) *DBLX) /DBLX 
S +U(I,J+1) 
UXY=(UX2-UX1) * (TR(ID, 4) -(NOY+J-1) *DBLY) /DBLY+UX1 
C..DEFINITION OF U ? 
TR(ID, 2) =UXY+TO 
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ELSE 
C..PTS OUTSIDE OF BLOCK 

DO 200 ID1=1,NOD1 

IF (ID.NE.ID1) THEN 

TSO=TR(ID1,2) 

KTS0=0 

DO 212 K=1,MNTSO+1 

IF (TSO.LE.TO-(K-1) *DTSO.AND.TSO.GT.TO-K*DTSO) KTSO=K 

212 CONTINUE 

C...TL=(ITL-1)*DTL, RL=(IRL-1) *DRL 

TL=TIME-TR(ID1,1) 

RL=SQRT((TR(ID,3)-TR(ID1,3) ) **2 


S +(TR(ID,4)-TR(ID1,4) ) **2) 
IF (RL. LE. (MNRL-1) *DRL.AND. 
$ TL. LE.TEND(KTSO) .AND.TL. LE. TEND (KTS0+1) ) THEN 


C..CALL MDU 
CALL MDU(UDATA,TSO,TL, RL, UMD) 
UIJ=UMD 
ELSE 
C..CALL FUNCTION 
IF (TL.LT.TEVAP(KTSO) ) THEN 
UIJ=UFEV(TSO,VO,BETA,TL,RL, KTSO, TEVAP, TEND) 
ELSE 
UIJ=UFRE (TSO, VO, BETA, TL,RL, KTSO, TEVAP, TEND) 
ENDIF 
ENDIF 
UOUT=UIJ+UO0UT 
ENDIF 
200 CONTINUE 
C..DEFINITION OF U=T-TO 
TR(ID,2)=UOUT+TO 
ENDIF 
C..NEW NOD 
NOD=NOD+1 
ELSE 
IF(TR(ID,1).GE.TIME+DTIME)GO TO 300 
ENDIF 
100 CONTINUE 
300 CONTINUE 


RETURN 
END 
Cw ec ee SUB. FOR FINDING A TEMP. (DIFFERENCE BETWEEN REAL 
Cesics SOLID SURFACE TEMP. AND INITIAL SOLID SURFACE TEMP.) 
Se eae s FUNCTION VARIABLES ARE: 
Cc. TSO: INITIAL SOLID SURFACE TEMP. 
es VO: INITIAL DROPLET VOLUME 
C. BETA: RATIO BETWEEN THE RADIUS OF THE WET AREA 
Cc. AND RADIUS OF THE EQUIVALENT SPHERICAL VOLUME 
C. 3 ig Fh LOCAL TIME IN WHICH TL=0 REPRESENTS THE 
oe MOMENT THE DROPLET STARTS TO EVAPORATE 
Cc. RL: LOCAL RADIUS WITH ORIGIN DEFINED AT 
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ANAANANAAN 


eeee#ee 


THE CENTER OF THE DROPLET 
SUBROUTINE MDU(UDATA,TSO,TL, RL, UMD) 
DIMENSION UDATA(MNTSO,MNVO,MNBETA,MNRL, MNTL) 
COMMON /DDATA/DTSO,DTL,DRL,TO,DTIME 
COMMON /NDATA/MNTSO,MNVO,MNBETA,MNRL,MNTL 
COMMON /NUMR/NRL,NOD,NDROP 
COMMON /GEOM/BLX, BLY, OBX,OBY,DBLX, DBLY 
COMMON /XY/NIX,NIY,NOX,NOY 

MNTSO=MN1 
MNVO=MN2 
MNBETA=MN3 
MNRL=MN4 
MNTL=MN5 


FIND ITSO TSO=TO-(ITSO-1) *DTSO 

DO 10 I=1,MNTSO+1 

IF(TSO.LE.TO-(I-1) *DTSO.AND.TSO.GT.TO-I*DTSO) THEN 
ITSO=I 

ENDIF 

CONTINUE 

FIND IVO 

IVO=1 

FIND IBETA 

IBETA=1 

FIND ITL TL=(ITL-1) *DTL 

ITL=INT (TL/DTL) +1 

FIND IRL RL=(IRL-1) *DRL 

DO 20 I=1,MNRL 

IF(RL.GE. (I-1) *DRL.AND.RL.LT.I*DRL) THEN 
IRL=I 

ENDIF 

CONTINUE 


shel ents INTERPOLATION IN T 


UT1=UDATA (ITSO,1,1,IRL,ITL) 
UT2=UDATA(ITSO,1,1, IRL, ITL+1) 
IF(UT1.GT.0..OR.UT2.GT.0.)PRINT *,ITSO,IRL,ITL 
UR1=(UT2-UT1) * (TL-(ITL-1) *DTL) /DTL+UT1 
UT1=UDATA (ITSO,1,1,IRL+1,ITL) 
UT2=UDATA(ITSO,1,1,IRL+1, ITL+1) 
IF(UT1.GT.0..OR.UT2.GT.0.)PRINT *,ITSO,IRL,ITL 
UR2=(UT2-UT1) * (TL-(ITL-1) *DTL) /DTL+UT1 
UT1=UDATA (ITSO0+1,1,1,IRL, ITL) 

UT2=UDATA (ITSO+1,1,1,IRL, ITL+1) 
IF(UT1.GT.0..OR.UT2.GT.0.)PRINT *,ITSO,IRL,ITL 
UR3=(UT2-UT1) * (TL-(ITL-1) *DTL) /DTL+UT1 
UT1=UDATA (ITSO0+1,1,1,IRL+1,ITL) 

UT2=UDATA (ITS0+1,1,1,IRL+1, ITIt1) 
IF(UT1.GT.0..OR.UT2.GT.0.)PRINT *,ITSO,IRL,ITL 
UR4=(UT2-UT1) * (TL-(ITL-1) *DTL) /DTL+UT1 


see tees INTERPOLATION IN R 


UMD1=(UR2-UR1) * (RL-(IRL-1) *DRL) /DRL+UR1 
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UMD2=(UR4-UR3) * (RL- (IRL-1) *DRL) /DRL+UR3 
IF(UMD1.GT.0..OR.UMD2.GT.0.)PRINT *,'UMD1,UMD2',UMD1,UMD2 
ae ek INTERPOLATION IN TSO 
UMD=- (UMD2-UMD1) * (TSO-(TO-(ITSO-1) *DTSO) ) /DTSO+UMD1 
IF(UMD.GT.0.)PRINT *,'UMD: TSO,RL,TL,UMD',TSO,RL,TL, UMD 
RETURN 
END 
Crewe at FUNCTION OF U FOR RL BEYOND (MNRL-1) *DRL DURING EVAPORATION 
FUNCTION UFEV(TSO,VO,BETA,TL,RL, KTSO, TEVAP, TEND) 
DIMENSION TEND (MNTSO) , TEVAP(MNTSO) 
COMMON /NDATA/MNTSO,MNVO,MNBETA, MNRL, MNTL 
PI=4.*ATAN(1.) 
Gre: se MACOR: 
SK=1.297E-3 
SR=2520.E-9 
SC=888.9 
C....DIFFUSIVITY OF MACOR 
ALPHAS=SK/ (SR*SC) 
C....LATENT HEAT (J/KG K) 
ALAMDA=2500.*1000. 
C....WATER DENSITY (KG/MM‘3) 
RHOW=997.6E-9 
C....LOCAL TIME IN FUNCTION (+1 TO AVOID THE TLO TO BE 0) 
TLO=TL : 
C....EVAP. TIME TEVAP1 
TEVAP1=TEVAP (KTSO) 
IF (TLO.LE.0.) THEN 
TL1=1. 
ELSE 
TL1=TLO 
ENDIF 
TSUB=TEVAP1-TL1 
UFEV=-2 . *V0* (RHOW/SR) * (ALAMDA/SC) / (8. * (PI*ALPHAS*TSUB) **1.5) 
$ *EXP(-RL**2/(4.*ALPHAS*TSUB) ) *COS (PI/2.*(1.-TLO/TEVAP1) ) 


RETURN 

END 
ele is see = FUNCTION OF U FOR RL BEYOND (MNRL-1)*DRL DURING THE 
lajste eta es RECOVERING PROCESS AND FOR TL EXCEEDING MNTL*DTL 


FUNCTION UFRE(TSO,V0,BETA,TL,RL,KTSO, TEVAP, TEND) 
DIMENSION TEND(MNTSO) , TEVAP (MNTSO) 

COMMON /NDATA/MNTSO,MNVO,MNBETA,MNRL, MNTL 
PI=4.*ATAN(1.) 


SK=1.297E-3 
SR=2520.E-9 
SC=888.9 
C....DIFFUSIVITY OF MACOR 
ALPHAS=SK/ (SR*SC) 
C....LATENT HEAT (J/KG kK) 
ALAMDA=2500.*1000. 
C....WATER DENSITY (KG/MM*3) 


RHOW=997.6E-9 


C....LOCAL TIME IN FUNCTION (+1 TO AVOID THE TLO TO BE 0) 
TLO=TL-TEVAP (KTSO) 
C...UFRE=? 


IF(TLO.LE.0.)TLO=.01 
UFRE=-2. *VO* (RHOW/SR) * (ALAMDA/SC) / (8. * (PI*ALPHAS*TLO) **1.5) 


$ *#EXP(-RL**2/(4.*ALPHAS*TLO) ) 
RETURN 


END — 
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Gia as segue 


RECOVERY PROCESS 
CONNECTED WITH TEST.FOR 
PARAMETER (JM=11,JE=4* (JM-1) ,NIMSL=JM+JE, KNT=5000) 
DIMENSION TN1(JM+JE) , THALF (JM+JE) 
DIMENSION TX(JM+JE) ,TR(JM+JE) , TB(JIM+JE,IM+JE) 
DIMENSION W(JM+JE-1,JM+JE-1) , FTIM(JM+JE-1,KNT) , F (JM+JE-1) 
DIMENSION W1(JM+JE-1,JM+JE-1, KNT) 
DIMENSION RHS(JM+JE-1) ,WO(JM+JE-1,IM+JE-1) 
DIMENSION RSO(JM+JE-1) ,WF(JM+JE-1) , TBM(JM+JE) 
DIMENSION PTT(50) 
DIMENSION IWK(NIMSL) , WK (NIMSL+NIMSL* (NIMSL-1) /2) 
DIMENSION TMD(300) 

COMMON /WORKSP/RWKSP 

COMMON /CONST/ ALFA, PI,KT,KKNT 

COMMON /GRID/ DRO,DTO,NDP,NP 

COMMON /INTG/ ATO,ERRREL, IRULE 

COMMON /SOLID/ SK,SC,SR 

COMMON /TSTEP/ NT 
CALL IWKIN(NIMSL+NIMSL* (NIMSL-1) /2) 


C...GET ALL INFORMATION FROM THE EVAPORATION PROCESS 


Cea e -AS 


10 


dee 


Cree 


THE INITIAL CONDITION OF THE RECOVERY PROCESS 
READ (15, *) ETIME, DTIME, TMPAIR, TSLDO, RR, RHK, KT, MMD 

READ (15, *)VO,BETA, IDMAT, DTO, NDP, NP, ATO, ERRREL, IRULE, ALFA, PI 
DO 10 I=1,JM+JE-1 

READ(15,*) (WO(I,J) ,J=1,IM+JE-1) 

CONTINUE 

READ(15,*) (TN1(I) , I=1,JM+JE) 
READ(15,*) (RSO(I) , I=1,JM+JE-1) 
READ(15,*) (RHS(I) , I1=1,JM+JE-1) 

DO 11 K=1,kKT 

READ(15,*) (FTIM(I,K) , 1=1,JM+JE-1) 

CONTINUE 

PRINT *,'NUMBER OF PLOTS?! 

READ(5,*)NPT 

PRINT *,'AT WHAT TIME (SECOND) '! 
READ(5,*) (PTT (NN) ,NN=1,NPT) 


TIME=-DTIME 
NT=0 

KKNT=1 
DTEMP=.005 


Conusone Op MATRIX 


90 
¢ Ue 


DO 90 I=1,JM+JE 
DO 90 J=1,JM+JE 
TB(I,J)=0. 
CONTINUE 


DO 150 J=2,JM+JE 
DO 155 JJ=2,JMt+JE 


TB(J,JJ)=RHK*WO(J-1,JJ-1) 
155 CONTINUE 
TB(J,J)=1.+RHK*WO (J-1,J-1) 
150 CONTINUE 
C...FOR J=1. B.C. T(I,1)=T(I,2) 
J=1 
TB(J,1)=1. 
TB(J,2)=-1. 
TR(J)=0. 
eWay) NORMALIZING THE MATRIX 
DO 133 I=1,JM+JE 
TBM (I)=TB(I,I) 
DO 133 J=1,JMtJE 
TB(I,J)=TB(I,J) /TBM(I) 
shake) CONTINUE 


CALL SOLVE1(JM+JE,TB,TR, TX, IWK, WK) 


TIME=TIME+DTIME 
Cinciemeri es: so a 0 bial) alGlele BIG)E WA STOVGATLTS « Glaih e ane! . KN t es RReerore cee <a 
TN10=TN1 (2) ; 
C....FOR PTS ON SOLID EXPOSED SURFACE 
DO 151 J=2,IM+JE 
SWJ=0. 
DO 157 JJ=2,IM+JE 
SWJ=SWJ+W0(J-1,J3-1) 
157 CONTINUE 
TR (J) =RHS (J-1) +(TSLDO-TMPAIR) * (1.+RHK*SWJ) 
151 CONTINUE 
Cri eae NORMALIZING THE MATRIX 
DO 143 I=1,JM+JE 
TR(I)=TR(I) /TBM(I) 
143 CONTINUE 


Ct eoeeeeeerere?€e? eoeeoeeeneee0e4eeereecskreesvseensv3eees3#*#ee3s5oe#eeeeege#eete ee @ @eeeee'e @ ese 019 * 2 


CALL SOLVE2 (JM+JE,TB,TR, TX) 


DO 300 J=1,IM+JE 
THALF (J) =.5*(TN1(J)+TX(J)) 
TN1(J)=TX(J) 
300 CONTINUE 
Clee UPDATE THE FLUX MATRIX 
KKT=KT+NT 
DO 417 KK=KKT,2,-1 
DO 417 II=1,JM+JE-1 
FTIM(II,KK)=FTIM(II,KK-1) 
417 CONTINUE 
DO 418 JJ=2,IM+JE 
Cc. cELUSES 
FTIM(JJ-1,1) =-RHK* (TN1 (JJ) +TMPAIR-TSLDO) 
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418 CONTINUE 
Cees eoeeoeeeee#segeeoeesqse#eee#eoe#eeee#ee#e#ese#esg$e#F#eee#e#e##28t##2eee#eeeeee# 
CALL NEXTR(NT,KKT,FTIM,RHS,TIME,W,F,WF,RSO,W1) 
C...TEST FOR CENTER POINT TEMPERATURE DIFFERENCE 
IF (NT.GT.1.AND.ABS(TN10-TN1(2)).LE.DTEMP)GO TO 3 
GO TO 1 


| PRINT *, 'THE SOLID Is RECOVERED AT TIME=',TIME 
STOP 
END 


SUBROUTINE NEXTR(NT,KT,FTIM,R,TIME,W,F,WF,RSO,W1) 
DIMENSION W(NP,NP) ,F(NP) , FTIM(NP, KT) 
DIMENSION WF(NP) ,R(NP) ,RSO(NP) 
DIMENSION W1(NP,NP,500) 
COMMON /CONST/ ALFA, PI,KNT,KKNT 
COMMON /GRID/ DRO,DTO,NDP,NP 
COMMON /INTG/ ATO,ERRREL, IRULE 
COMMON /MATER/ LMAT 
DO 149 II=1,NP 
R(II)=0. 
149 CONTINUE 
KKT=KT 
DO 150 IT=1,KKT 
C SELECT PROPER WEIGHTS 
TO=IT*DTO 
azn 

IF (IT.LE.KKNT) THEN 

DO 151 I=1,NP 

DO 151 J=1,NP 

W(I,J)=W1(I,J,IT) 

151 CONTINUE 
ELSE 


CALL WEIGHT (TO,W,RSO) 


DO 152 I=1,NP 
DO 152 J=1,NP 
W1(I,J,IT)=W(I,J) 
152 CONTINUE 
ENDIF 
C PICK UP CORRESPONDING FLUX 
DO 100 I=1,NP 
F(I)=FTIM(I,IT) 
100 CONTINUE 
C SUMMATION IN SPACE 
CALL MULT(NP,NP,W,F,WF) 
C SUMMATION IN TIME 
CALL ADD(NP,WF,R) 
150 CONTINUE 
KKNT=KKT 
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20 


50 


Cc 
(@: 


RETURN 
END 

SUBROUTINE MULT(M,N,W,F,WF) 
REAL W(M,N),F(N) ,WF(N) 

DO 50 I=1,M 

WF(I)=0.0 

DO 50 K=1,N 

WF (I) =WF(I)+W(1I,K) *F(K) 
CONTINUE 

RETURN 

END 


SUBROUTINE ADD(NP, WF,R) 
REAL WF(NP) ,R(NP) 

DO 50 I=1,NP 
R(I)=R(1I)+WF(I) 
CONTINUE 

RETURN 

END 


SUBROUTINE WEIGHT (TO,W,RSO) 
REAL W(NP,NP) ,RSO(NP) 

REAL QDAG, BSIOE, ERF 
EXTERNAL Fl 


COMMON /CONST/ ALFA, PI,KNT, KKNT 
COMMON /GRID/ DRO,DTO,NDP,NP 
COMMON /INTG/ ATO,ERRREL, IRULE 
COMMON /RRO/ R,RO 

COMMON /MATER/ LMAT 


DO 50 I=1,NP 
DO 50 J=1,NP 

IF(J.EQ.1) THEN 
DRO=2.*RSO (J) 
ELSE 

IF(J.EQ.NP) THEN 

DRO=DRO 
ELSE 

DRO=.5* (RSO(J+1) -RSO(J-1) ) 
ENDIF 7 
ENDIF 

=RS0 (I) 
RO=RS0 (J) 

PRINT *,DRO,R,RO 


SETAU PaGioeenOk o2 0 < DITO 


ERRABS=0.0 
IF (TO.LT.DTO) THEN 


W(I,J)=QDAG(F1,ATO,.5,ERRABS , ERRREL, IRULE, RESULT, ERREST) 


W(I,J)=RESULT 
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END IF 
C SET UP (W) FOR TO > DTO 

IF(TO.GE.DTO) THEN 

ARG1=DRO/SQRT (16*ALFA*TO) 

ARG3=(R-RO) **2/ (4*ALFA*TO) 

ARG4=0.5*R*RO/ (ALFA*TO) 

IF(I.NE.J) THEN 

W(I,7)=(1.0/SQRT(4.0*PI*ALFA) ) *RO* (TO**(-1.5))* 

& BSIOE(ARG4) *EXP(-ARG3) *DRO*DTO 

ELSE 

W(1I,J)=(R/TO) *BSIOE(ARG4) *ERF (ARG1) *DTO 

END IF 

END IF 


50 CONTINUE 
RETURN 
END - 


FUNCTION F1(TO) 
REAL BSIOE,ERF 

COMMON /CONST/ ALFA, PI,KNT,KKNT 
COMMON /GRID/ DRO,DTO,NP 
COMMON /RRO/ R,RO 


ARG1=DRO/SQRT (16*ALFA*TO) 
ARG3=(R-RO) **2/ (4*ALFA*TO) 
ARG4=0.5*R*RO/ (ALFA*TO) 
E=0.0 
IF (ARG3.LE.100.0) E=EXP(-ARG3) 
IF(R.NE.RO) THEN 
F1=(1.0/SQRT(4.0*PI*ALFA) ) *RO* (TO**(-1.5))* 
& BSIOE(ARG4) *E*DRO 

ELSE 
F1=(R/TO) *BSIOE(ARG4) *ERF (ARG1) 
END IF 
RETURN 
END 

SUBROUTINE SOLVE1(N,A,B,X, IWK, WK) 


eLey 4. SUBROUTINE FOR SOLVING MATRIX 
C... 'LINK®XXX, SOLVEPIMSL/LIB' 97: 
DIMENSION A(N,N),B(N) , X(N) 
DIMENSION IWK(N) ,WK(N+N* (N-1) /2) 
INTEGER N 
COMMON /WORKSP/RWKSP 
CALL L2NRG(N,A,N,A,N,WK, IWK) 
RETURN 
END 
ae 
SUBROUTINE SOLVE2(N,A,B,X) 
C.....SUBROUTINE FOR SOLVING MATRIX 
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C...! LINKEXXRISOLVE, IMSL/LIB"..... 
DIMENSION A(N,N) ,B(N) , X(N) 
INTEGER N 
DO 100 J=1,N 
X(J)=0.0 
DO 100 I=1,N 
X(J)=X(J)+A(J,1I) *B(I) 

100 CONTINUE 
RETURN 
END 
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C....-PROGRAM FOR SOLVING A CONDUCTION HEAT TRANSFER PROBLEM IN A 


C....LIQUID DROPLET LAID ON THE SURFACE OF AN INFINITE LARGE SOLID. 
‘schog eee THE BOUNDARY ELEMENT METHOD IS USED FOR THE SOLID 
Gente! AND THE INTEGRAL CONTROL VOLUME METHOD IS USED FOR THE 
O3..... <xeLIQUID 
ase eae CRANK-NICHOLSON SCHEME IS USED 
Os 5 se IM: NUMBER OF POINTS IN LIQUID DROP IN VERTICAL DIRECTION 
Cees JM: NUMBER OF POINTS IN LIQUID DROP IN HORIZATAL DIRECTION 
Gs... re. JE: NUMBER OF POINTS ON EXPOSED SOLID SURFACE 
eo ee KNTO: NUMBER OF RECOLLECTION TIME STEPS 
.... Dr NIMSL: TOTAL GRID POINTS 
PARAMETER (IM=16,JM=11,JS=JM-1,JE=4*JS, 
$ KNTO=60, NIMSL=IM*JM+JE) 
DIMENSION TN1(IM*IM+JE) , TR(IM*IM+JE) , THALF (IM*IM+JE) 
DIMENSION TX(IM*JM+JE) , TA(IM*IM,4) , TB(IM*IM+JE, IM*JM+JE) 
DIMENSION W(JM+JE-1,JM+JE-1) , FTIM(JM+JE-1,KNTO) , F(JM+JE-1) 
DIMENSION W1(JM+JE-1,IM+JE-1,KNTO) 
DIMENSION RHS(JM+JE-1) ,CT(JM) ,WO(JM+JE-1,JM+JE-1) 
DIMENSION EWIS(JM) ,SX(JM) ,SX2 (JM) 
DIMENSION R(IM+1,JM+1) ,Z(IM+1,JM+1) ,DZ(JM) ,RSO(5*JS) 
| DIMENSION RO3(IM) ,Z03(IM) , PTT(50) 
DIMENSION IWK(NIMSL) , WK(NIMSL+NIMSL* (NIMSL-1) /2) 
DIMENSION TMD(300) , TPLOT( (IM-1) * (JM-1) +1) 
DIMENSION TBS (JM+JE,JM+JE) , TRS (JM+JE) , TXS (JM+JE) 
COMMON /SUBS/TIME, DTIME, TMPAIR, TSLDO, RR, KNT, RHK 
COMMON /WORKSP/RWKSP 
COMMON /CONST/ ALFA, PI 
COMMON /GRID/ DRO,DTO,NDP,NP 
COMMON /INTG/ ATO,ERRREL, IRULE 
COMMON /SOLID/ SK,SC,SR 
COMMON /WATER/ COEFKW, CPW, RHOW 
COMMON /TSTEP/ NT 
Can 
PI=4.*ATAN(1.0) 
TRULE=2) 
ATO=0.001 
ERRREL=0.001 
DTIME=.1 
NTIME=500 
C....CHOICE OF THE MATERIAL FOR THE SOLID SURFACE 
PRINT *,' SOLID = ? (TEST=0, ALUMINUM=1, MACOR=2, 
$  STEEL=3,MARBLE=4,BAKELITE=5) ' 
. READ(5,*) LMAT 
C....INITIAL SETTING THE DROPLET SIZE BY KNOWING VO AND BETA 
PRINT *,' INITIAL VOLUME OF THE DROPLET= ?!' 
READ(5,*) VO 
PRINT *,' BETA= ?! 
READ(5,*) BETA 
C...RR: RADIUS OF THE WETTED AREA 
C...AA: APEX OF THE DROPLET 


RR=BETA*EXP(1./3.*LOG(.75*VO/PI) ) 
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UU=EXP(1./3.*LOG(3.*VO/PI+SORT (RR**6+9. /PI**2*VO*#2) )) 
AA=UU-RR**2/UU 


ay 
MMD=INT (10*RR) 
en 
VO0=VO 
AAO= 


C. .).'.' VERIFY 
VO=PI*AA/6.%*(3.*RR**2+AA¥**2) 


PRINT *,'RR=',RR,' AA=',AA,' VO=',VO,!' DT=',DTIME 
C...AIR TEMP. < 

TMPAIR=20. 
C...INITIAL TEMPERATURE OF THE LIQUID DROP 

TEMPO=20. 


C...INITIAL TEMPERATURE OF THE SOLID 
PRINT *,'INITIAL SOLID TEMPERATURE: ' 
READ (5, *) TSLDO 
C...ALL TEMPERATURES USED HERE, TN1 AND TX, ARE DELTAT : 
U=TEMP-TMPAIR ! 
TAIR=TMPAIR-TMPAIR 
C...SOLID CONDUCTIVITY (SK:W/MM K; SR:KG/MM*3; SC:J/KG k) 
IF (LMAT.EQ.0) THEN 
PRINT *,'SK(W/MM K), SR(KG/MM‘*3), SC(J/KG K)' 
READ (5,*)SK,SR,SC 
ENDIF 
Cr eee . ALUMINUM: 
IF (LMAT.EQ.1) THEN 
SK=180.E-3 
SR=2771.E-9 


IF (LMAT.EQ.2) THEN 
SK=1.297E-3 
SR=2520.E-9 


IF (LMAT.EQ.3) THEN 
SK=12.E=-3 
SR=8100.E-9 


IF (LMAT.EQ.4) THEN 
SK=2.8E-3 
SR=2600.E-9 


IF (LMAT.EQ.5) THEN 
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SK=.23E-3 
SR=1300E-9 
SC=1600. 
ENDIF 


C...CONVECTIVE HEAT TRANSFER COEFF. FOR LIQUID SURFACE (W/MM“%2 Kk) 


COEFFH=(.098*TSLDO0+2.79) *1.E-6 
C...HEAT TRANSFER COEFF. FOR EXPOSED SOLID SURFACE 
SH=COEFFH 
C...RATIO OF SOLID H AND K : H/K 
RHK=SH/SK 
C...DIFFUSIVITY OF SOLID (MM*2/S) 
ALFA=SK/ (SR*SC) 


C...WATER DENSITY (KG/MM‘*3) 
RHOW=997.6E-9 

C...-WATER CONDUCTIVITY (W/MM K) 
COEFKW=.613E-3 

C...WATER SPECIFIC HEAT(J/KG K) 


CPW=4179. 
C...DIFFUSIVITY OF WATER: ALPHA=KW/(RHO*CP) (MM‘2/S) 
ALPHA=COEFKW/ (RHOW*CPW) 
C...RATIO BETWEEN LIQUID AND SOLID CONDUCTIVITES 
RKLS=COEFKW/SK 


C...AIR SPECIFIC HEAT (J/KG K) 
CA=1009. 
C...AIR PRESSURE (BAR) 
PA=1.014 
C..XA IS DETERMINED WHEN PJ IS EVALUATED AT AIR TEMP. 
« ' - XA=EXP(11.7789-3876.0717/ (TMPAIR+229.4483) )/PA 
XA=.012 
C....LAMDA (LATENT HEAT): (J/KG) 
Cc ALAMDA=(2509.-2.51*TJ) *1000. 


ALAMDA=2250.*1000. 
ALAMD=2500.*1000. 
ret cy 3 ..FOR PLOT INFORMATION 
PRINT *,'HOW MANY PLOTS YOU NEED?! 
READ (5,*)NPT 
PRINT *,'AT WHAT TIME (SECOND) ' 
READ(5,*) (PTT (NN) ,NN=1,NPT) 
C...INITIAL SETTING THE DROPLET TEMP. 
DO 5 I=1,IM 
DO 5 J=1,0M 
IJ=(I-1) *TM+I 
TN1 (IJ) =TEMPO-TMPAIR 
5 CONTINUE 
C....CONTACT TEMPERATURE 
CALL CTEMPO(IM,JM,JE,TMPAIR, TN1,TSLDO,CT,TTO) 
C....SETTING INITIAL TEMP. FOR FAKE PTS. (TN1(I=IM) ) 
DO 51 J=1,JIMtJE 
IJ=(IM-1) *IM+d 
TN1(1IJ)=TSLDO-TMPAIR 
ENDIF 
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52 


CONTINUE 


C...MAXIMUM TEMP. LIMIT 


TLIMIT=99.5-TMPAIR 
RRO=AA 


QO=RHOW*VO*ALAMD/ (PI*RR*¥*2) 


TEO=RRO*RHOW* CPW/COEFFH 
TFINT=0. 

TF1=0. 

TF2=0. 


C...-FOR AN INITIAL TEMP SETUP 


999 


KL=0 
KKNT=0 
TIME=-0.5*DTIME 
IF (KKNT.EQ. 0) THEN 
KNT=10 

ELSE 

KNT=KNTO 

AA=AAO 

VO=Vo0 

ENDIF 

NT=0 


DO 53 I=1,JM+JE-1 
RHS (I) =0. 
CONTINUE 


eeese#eee#sgreer#+r»reeeee eoeeeeee#eeetese 


DO 1 NT=1,NTIME 
NT=NT+1 


C...SET NEW GRID POINTS BASED ON NEW RR,AA 
CALL RZ90(IM,JM,RR,AA,R,Z,JS,RSO) 
CALL TRIO(IM,JM,R,Z,RO3,2Z03) 

CALL DZ0(IM,JM,R,Z,RO3,2Z03,DZ) 


81 


eoeee53s5oiteeeeo0ee3evete3e#enreoeese3ee#ee8e8eee8 @ 


IF (NT.EQ.1.AND.KKNT.EQ.0) THEN 


DTIME=.1 
TTT=DTIME*.99 
NDP=JM-1 
NP=JM+JE-1 
DTO=DTIME 
CALL WEIGHT (TTT,WO,RSO) 
ENDIF 
CONTINUE 


C...SET OFMATRIZ 


ee) 


DO 90 I=1,IM*JIM+JE 

DO 90 J=1,IM*JIM+JE 

TB(I,J)=0. 
CONTINUE 


C...-FOR LIQUID: INNER POINTS 


DO LOOTAL=27, IM=<1 
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DO 100 J=2,JM-1 
CALL NEWSO(IM,JM,R,Z,RO3,2Z03, 

$ I,J,RN,2ZN,RE, ZE, RW, ZW,RS,ZS,RO, ZO) 
ALAB=SQRT ( (R(I+1,J+1)-R(I,J+1) ) **2+(Z(I+1,J+1) -Z(I,J+1) ) **2) 
ALBC=SQRT((R(I,J+1)-R(1I,J) ) **2+(Z(1I,J+1)-Z(1I,J) ) **2) 
ALCD=SQRT((R(I,J)-R(I+1,7) ) **2+(Z(I,J7)-Z(I+1,7) ) **2) 
ALDA=SQRT ( (R(I+1,J7) -R(I+1,J+1) ) **2+(Z(I+1,7) -Z(I+1,J+1) ) **2) 
AREA=.25* (ALAB+ALCD) * (ALBC+ALDA) 
DXN=SQRT ( (RN-RO) **2+(ZN—-ZO) **2) 
DXW=-SQRT ( (RO-RW) **2+ (ZO-ZW) **2) 
DXS=-SQRT ( (RO-RS) **2+(ZO-ZS) **2) 

DXE=SQRT ( (RE-RO) **2+(ZE-ZO) **2) 
IJ=(I-1) *IM+I 
Orr. 
CONST=.25*DTIME*ALPHA/ (RO*AREA) 
AIM1J=CONST*ALBC# (R(1I,J+1)+R(1I,J) ) /DXN 
AIJP1=CONST*ALAB* (R(I+1,J+1)+R(I,J+1) ) /DXE 
AIJM1=-CONST*ALCD* (R(1I,J)+R(I+1,J) ) /DXW 
AIP1J=-CONST*ALDA* (R(I+1,J)+R(I+1,J+1) ) /DXS 
AIJ=- (AIM1J+AIJM1+AIJP1+AIP1J) 
Cee: 
TB(IJ,IJ+1)=-AIJP1 
TB(IJ,IJ+JM) =-AIP1J 
TB(IJ,IJ-JM) =-AIM1J 
TB(IJ,IJ-1)=-AIJM1 
5 Ga faa 09 Oo ye 
Coe 
IF (NT.EQ.1) THEN 
PAG: th) =A bie): 
TA (Id, 2SALRIO 
TA (IJ,3)=AIM1J 
TA(IJ,4)=AIJM1 
ENDIF 
TR(IJ)=TA(IJ,1) *TN1(IJ+1)+TA(IJ,2) *TN1 (IJ+dM) 
$+TA(IJ,3) *TN1(IJ-JM)+TA(IJ,4) *TN1(IJ-1) 
$+(1.-(TA(IJ,3)+TA(IJ,4)+TA(IJ,1)+TA(IJ,2))) *TN1(IJ) 
C..SAVE FOR NEXT TIME STEP 
TA(IJ;1)=AIJP1 
TA(IJ,2)=AIP1I 
TA(IJ,3)=AIM1J 
TA(IJ,4)=AIJM1 
(eve 
100 CONTINUE 
C....FOR J=JM: EDGE POINTS, TRIANGULAR ELEMENTS 
J=IM 
DO 110 I=2,IM-1 
CALL NEWSO(IM,JM,R,Z,RO3,2Z03, 

S I,J,RN,ZN,RE, ZE, RW, ZW, RS, ZS,RO, ZO) 
ALCA=SQRT((R(I,J+1) -R(I+1,d) ) **2+(Z(I,J+1) —Z (I+1,J) ) **2) 
ALAB=SOQRT((R(I,J+1)-R(I,J) ) **2+(Z(1I,J+1) -Z (I,J) ) **2) 
ALBC=SQRT((R(I+1,J)-R(I,d) ) **2+(Z(I+1,J) -Z(I,J) ) **2) 
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$ 
$ 


AREA=.25* (ALCA+ALAB) * (ALBC+ALAA) 
DXN=SQRT ( (RN-RO) **2+(ZN-ZO) **2) 
DXW=-SQRT( (RO-RW) **2+(ZO-ZW) **2) 
DXS=-SORT( (RO-RS) **2+(ZO-ZS) **2) 
IJ=(I-1) *OM+T 


CONST=.25*DTIME*ALPHA/ (RO*AREA) 

AIJP1=0. 

AIP1J=-CONST*ALCA* (R(I+1,J)+R(1I,J+1) ) /DXS 
AIM1J=CONST*ALAB* (R(I,J+1)+R(I,J) ) /DXN 
AIJM1=-CONST*ALBC* (R(I,J)+R(I+1,J) ) /DXW 
AIJ=-(AIM1J+AIJM1+AIJP1+AIP1J) 


TB(IJ,IJ+1)=-AIJP1 
TB(IJ,IJ+JIM) =-AIP1T 
TB(IJ,IJ-JM) =-AIM1I 
TB(IJ,IJ-1) =-AIJM1 
TB(IJ,IJ)=1.-AIT 
IF(NT.EQ.1) THEN 
TA(IJ,1)=AIJP1 
TA(IJ,2)=AIP1I 
TA(IJ,3)=AIM1T 
TA(IJ,4)=AIJM1 
ENDIF 
TR(IJ)=TA(IJ,2) *TN1(IJ+JM) 


+TA(IJ,3)*TN1(IJ-JM)+TA(IJ,4) *TN1(IJ-1) 


+(1.-(TA(IJ,3)+TA(IJ,4)+TA(IJ,1)+TA(IJ,2)))*TN1(IJ) 


C..SAVE FOR NEXT TIME STEP 


Cars 


TA(IJ,1)=AIJP1 
TA(IJ,2)=AIP1I 
TA(IJ,3)=AIM1T 
TA(IJ,4)=AIJM1 


CONTINUE 


C..FOR I=1:) LIQUID —- AIR INTERFACE 


Crue 


I=1 
DO 120 J=2,JUM 


- PROPERTIES BASED ON THE DROPLET SURFACE TEMP. TJ (DEGREE C) 


TJ=.5*(TN1(J)+TN1(J+JM) )+TMPAIR 


C....SURFACE PRESSURE (BAR). 


Cir 


-- 10<TJ<100, OTHERWISE SX<0. AND DVDT<O. 


IF(TJ.LE.10..AND.TJ.GE.100.) PRINT *,'TJ IS OUT RANGE FOR SX! 
IF(TJ.GE.99.5) THEN 

PJ=1. 

ELSE 

PJ=EXP (11.7789-3876.0717/ (TI+229.4483) ) 

ENDIF 

XJ=PI/PA 

SX (J) =(XJ-XA) /(1.-XJ) 


C...LEWIS NUMBER 
EWIS (J) =6.0* (-.0005279*TJ+1.152) 


IJ=J 
IF(J.EQ.JM) THEN 
RM=RO3 (I) 
ZM=Z03 (I) 
=RO3 (I+1) 
ZP=Z03 (I+1) 
ELSE 
RM=.25*(R(I,J)+R(I,J+1)+R(I+1,J+1)+R(I+1,J) ) 
ZM=. 25% (Z(1,3)+Z(1,0+1) +2 (11, 0+1)42(I+1,J),) 
RP=.25*(R(I+1,J)+R(I+1,J3+1)+R(I+2,J3+1)+R(I+2,J)) 
ZP=.25*(Z(I+1,J)+Z(I+1,J+1)+Z(I+2,JI+1)+Z(I+2,T)) 
ENDIF 
DZ1=SQRT ( (RM-RP) **2+(ZM-ZP) **2) 
TB(IJ,IJ+JM) =1.-.5*COEFFH*DZ1/COEFKW 
TB(IJ,IJ)=-1.-. 5*COEFFH*DZ1/COEFKW 
C...CONST1 IS FROM THE B.C. 
ALAMDA=(2509.-2.51*TJ) *1000. 
CONST1=.624*DZ1*COEFFH*ALAMDA*EWIS (J) *SX(J) / (COEFKW*CA) 
CNL 
TR (IJ) =CONST1-TAIR*DZ1*COEFFH/COEFKW 
120 CONTINUE 
C..FOR I=IM: LIQUID - SOLID INTERFACE 
C...COUPLED BOUNDARY ELEMENT AND CONTROL VOLUME METHODS 
I=IM 
PRINT *,'TIME=',TIME 
en ee ee ad ital ee ee oe SBA RA a iy: 
DO 130 J=2,JM 
IJ=(I-1) *IM+I 
DO 135 JJ=2,IM 
IJ1=(I-1) *IM+IT 
IJ2=(I-2) *IM+IT 
IF(IJ1.NE.IJ) THEN 
TB(IJ,1IJ1)=RKLS*W0 (J-1,JJ3-1) /DZ (JJ) 
TB(IJ,1IJ2)=-RKLS*WO0 (J-1,JJ3-1) /DZ (JT) 
ELSE 
TB(IJ,1J1)=.5+RKLS*WO (J-1,J-1) /DZ (JJ) 
TB(IJ,1J2)=.5-RKLS*WO (J-1,J-1) /DZ (JJ) 
ENDIF 
135 CONTINUE 
DO 136 JJ=JM+1,IM+JE 
IJ3=(I-1) *IM+IT 
TB (IJ, 133) =RHK*WO (J-1,J3-1) 
136 CONTINUE 
TR (IJ) =RHS (J-1) 
SWJ=0. 
DO 137 JJ=2,IM+JE 
SWJ=SWI+WO (J-1,J3-1) 
137 CONTINUE 
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TR (IJ) =TR(IJ)+(TSLDO-TMPAIR) * (1.+RHK*SWJ) 
130 CONTINUE 
C....FOR PTS ON SOLID EXPOSED SURFACE 
DO 150 J=JM+1,JM+JE r 
IJ=(I-1) *IM+T 
TB(IJ,1J)=1.+RHK*WO (J-1,J-1) 
DO 155 JJ=2,IM 
TJ1=(I-1) *IM+JT 
IJ2=(I-2) *IM+II 
TB(IJ,1IJ1)=RKLS*W0 (J-1,JJ-1) /DZ (JJ) 
TB(IJ,1J2)=-RKLS*WO(J-1,JI-1) /DZ (JJ) 
155 CONTINUE 
DO 156 JJ=JM+1,JM+JE 
IJ3=(I-1) *IM+II 
IF (IJ3.NE.IJ)TB(IJ,1IJ3) =RHK*WO(J-1,JJ-1)- 
156 CONTINUE 
TR(IJ)=RHS (J-1) 
SWJ=0. 
DO 157 JJ=2,IM+JE 
SWJ=SWI+WO0 (J-1,J3-1) 
157 CONTINUE | 
TR(IJ)=TR(IJ)+(TSLDO-TMPAIR) * (1.+RHK*SWJ) 
150 CONTINUE 


C...FOR J=1. B.C. T(I,1)=T(I,2), SYMMETRIC TO Z-AXIS 


J=1 
DO 140 I=1,IM 
IJ=(I-1) *OM+T 
TB(IJ,IJ)=1. 
TB(IJ,IJ+1)=-1. 
TR(IJ)=0. 

140 CONTINUE 

Chee NORMALIZING THE MATRIX 

DO 143 I=1,IM*JIM+JE 
TBM=TB (I,I) 
DO 144 J=1,IM*IM+JE 
TB(I,J)=TB(I,J) /TBM 


144 CONTINUE 
TR(I)=TR(1I) /TBM 
143 CONTINUE 


C...TEMPERATURE CALCULATION BY SOLVING TB*TX=TR 


C....SOLVE FOR EVERY STEP AFTER PRE-SETTING A TEMP. 


IF (NT.NE.1.OR.KKNT.NE.0O) THEN 


Cc coeee#eeee# eoeeeeseeeseee#egeeses# eeeeees7eeee3ee3ses#e##ee8e#e &@& @ &© @ @ eeeee#es: 


CALL SOLVE(IM*JM+JE,TB,TR, TX, IWK, WK) 
Kir) 
C....AVERAGE BETWEEN CURRENT AND PREVIOUS STEPS 

DO 300 IJ=1,IM*IM+JE 
THALF (IJ) =.5*(TN1(IJ)+TX (IJ) ) 
TN1 (IJ) =TX (IJ) 

300 CONTINUE 
JU=(IM-2) *JM 


DISTRIB 
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1 L 


310 


JL=(IM-1) *JM 
MDDATA DRI=1. 
TMD (1)=.25* (THALF (JU+1) +THALF (JU+2) +THALF (JL+1) 
+THALF (JL+2))+TMPAIR 

DO 910 I=1,MMD 
DO 911 J=2,IM+JE-1 
IF (I.GE.RSO(J-1) .AND.I.LT.RSO(J) ) THEN 
IF (I.LT.RR) THEN 
TMD1=.5* (THALF (JU+J) +THALF (JL+J) ) 
TMD2=. 5* (THALF (JU+J+1) +THALF (JI4J+1) ) 
TMD (I+1) =(TMD2-TMD1) * (I-RSO(J-1) )/ 
(RSO(J)-RSO(J-1) )+TMD1+TMPAIR 
ENDIF 
IF (I.EQ.RR) THEN 
TMD1=.5* (THALF (JU+JM) +THALF (JL+JM) ) 
TMD2=THALF (IM*JM+1) 
TMD (I+1) =(TMD2-TMD1) * (I-RSO(J-1) )/ 
(RSO(J)-RSO(J-1) )+TMD1+TMPAIR 
ENDIF 
IF (I.GT.RR) THEN 
TMD1=THALF ( (IM-1) *JM+J) 
TMD2=THALF ( (IM-1) *JM+J+1) 
TMD (I+1) =(TMD2-TMD1) * (I-RSO(J-1) )/ 
(RSO(J) -RSO(J-1) )+TMD1+TMPAIR 
ENDIF 
ENDIF 
CONTINUE | 
TEDGE=THALF ( IM*JM) + (THALF (IM*JM+1) 
-THALF (IM*JM) ) * (RR-RSO (JS) ) 

/ (RSO(JS+1) -RSO(JS) )+TMPAIR 
T3R=THALF (IM*JM+2.*JS)+TMPAIR 
IF (I.GE.RSO(JM+JE-1) ) THEN 
TD (I+1)=TSLDO 
ENDIF 
CONTINUE 
FORCING ALL TEMP. UNDER THE LIMIT 
DO 310 I=1,2 
DO 310 J=1,JM 
IJ=(I-1) *IM+d 
IF(TN1(IJ) .GE.TLIMIT) TN1(IJ)=TLIMIT 
CONTINUE 


C....CAL NEW DROPLET VOLUME AND AA 


CALL SURF (DTIME,IM,JM,RR,VO,AA,R,Z,DVDT, 
EWIS , COEFFH, SX, RHOW, CA) 


C... UPDATE THE FLUX MATRIX 


417 


DO 417 KK=KNT,2,-1 

DO. 417 II=1 ,JMtJE-1 
FTIM(II,KK) =FTIM(II, KK-1) 
CONTINUE 


eoeeee 
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C... COMPUTE THE FLUX AT CURRENT TIME: FTIM(II,1) 
DO 418 JJ=2,IM+JE 
J1=(IM-1) *IM+JIT 
J2=(IM-2) *IM+JT 
I=IM 
C... FLUXES 
IF (JJ.LE.JM) THEN 
FTIM(JJ-1,1)=RKLS* (TN1(J2)-TN1(J1) ) /DZ (JJ) 


S +RHK* (TSLDO-TMPAIR) 
ELSE 
FTIM(JJ-1, 1) =-RHK* (TN1 (J1) +TMPAIR-TSLDO) 
ENDIF 


418 CONTINUE 
C...FOR SOLID SURFACE POINTS 
CALL NEXTR(NT,KNT,FTIM,RHS, TIME,W,F,WF,RSO,W1) 
C...DETERMINE IF THE DROPLET IS COMPLETELY EVAPORATED 
IF (DVDT*DTIME.GT.VO.OR.AA.LT. .01) THEN 
DO 700 J=1,IM+JE 
IF(J.LE.JM) THEN 
TN1(J)=.5* (TN1( (IM-1) *JM+J)+TN1((IM-2) *JM+J) ) 
ELSE 
TN1(J)=TN1((IM-1) *JM+J) 
ENDIF 
700 CONTINUE 
ETIME=TIME 
PRINT *,'THE DROPLET IS GONE AT TIME=',ETIME 
ENDIF 
IF (NT.EQ.KNT+1.AND.KKNT.EQ.0) THEN 
DO 380 KK=KNT,1,-1 
DO 380 II=1,JM+JE-1 
FTIM(II,KK)=0. 
380 CONTINUE 
KKNT=1 
TIME=0. 
DTIME=.5 
DTO=DTIME 


CO etetale enc tetenene eoee7we#50jeees ee0ee7e0eeeekeeee8e28282028830 8080 @ @e@eoeoee8eo_ oe@eoeoeo#e#oe#et?ee¢ @ 6 8)8 8 se 


2 STOP 
END 
C...SUB. ASSIGNING THE CENTER PTS. FOR THOSE CELLS AT EDGE WHICH 
C...ARE TRIANGLES. 
SUBROUTINE TRIO(IM,JM,R,Z,RO3,2Z03) 
DIMENSION R(IM+1,JM+1) ,Z(IM+1,JM+1) 
DIMENSION RO3(IM) ,Z03(IM) 
J=IM 
DO 10 I=1,IM 
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C 


10 


$ 


A=(.5*(Z(I,J0)+Z(I,Jd+1) )-Z(I+1,J))/ 
(.5*(R(I,J)+R(I,J+1) )-R(I+1,J) ) 

B=(.5*(Z(I+1,J7)+Z(1I,Jd+1))-Z(1,J5))/ 
(.5*(R(I+1,7)+R(I,J+1))-R(I,J) ) 

IF (ABS (A-B) .GT.1.E-7) 

RO3 (I) =(Z(I,J)-Z(I+1,J7)+A*R(I+1,J)-B*R(I,J)) /(A-B) 

Z03 (I) =Z(1I,J)+B*(RO3(I)-R(I,J)) | 

CONTINUE 

RETURN 

END 


‘SUBROUTINE NEWSO(IM,JM,R,Z,RO3,Z03,1,J3,RN,ZN,RE,ZE, 


RW, ZW,RS,ZS,RO,ZO) 


, 


C....-SUBROUTINE FOR ASSIGNING THE POSITIONS AT NORTH, EAST, WEST, 
Cin SOU Lay OF (oA, CELL: 


C....THE LEFT UPPER CORNER IS THE POSITION R(I,J) AND Zan 


ehere ve he iene 


DIMENSION R(IM+1,JM+1) ,Z(IM+1,JM+1) 

DIMENSION RO3 (IM) , Z03 (IM) 

IF (J.EQ.JM) THEN 

RO=RO3 (I) 

ZO=Z03 (I) 

RN=RO3 (I-1) 

ZN=Z03 (I-1) 

RS=RO3 (I+1) 

ZS=Z03 (I+1) 
RW=.25*(R(I,J)+R(I,J-1)+R(I+1,J-1)+R(I+1,7)) 
ZW=.25*(Z(1,J0)+Z(1,J-1)+Z(I+1,J-1)+Z(I+1,7)) 
ELSE 
RO=.25*(R(1I,J)+R(1I,J+1)+R(I+1,J+1)+R(I+1,J)) 
ZO=,25*(Z(L,0)+Z(L,d+1)+Z(It+1,3+1)+Z(I+1,d)) 
RE=.25*(R(I,J+1)+R(1I,J+2)+R(I+1,J3+2)+R(I+1,J+1) ) 
ZE=.25*(Z(I,J+1)+Z(1,J+2)+Z(I+1,J+2)+Z(I+1,J+1) ) 
RN=.25*(R(1I,J+1)+R(I-1,J+1)+R(I-1,J7)+R(I,J)) 
ZN=.25*(Z(1,Jd+1)+Z(I-1,3+1)+Z(I-1,3)+Z(I,J)) 
RW=.25*(R(I,J)+R(1I,J-1)+R(I+1,J-1)+R(I+1,J)) 
ZW=.25*(Z(1,J)+Z(1,J-1)+Z(I+1,J-1)+Z(I+1,7) ) 
RS=.25*(R(I+1,J)+R(I+1,J+1)+R(I+2,J+1)+R(I+2,J)) 
ZS=.25*(Z(I+1,J0)+Z(I+1,J+1)+Z(I+2,JT+1)+Z(1+2,T)) 
ENDIF 


SUBROUTINE RZ90(1IM,JM,RR,AA,R,Z,JS,RSO) 


SUBROUTINE FOR ASSIGNING THE GRID POINTS..... 
R(I,J), 2(1,3): I=1,IM+1, J=1,IM+1 
IMAGINARY POINTS: 1) R(1,J), Z(1,3) 


2) R(IM+1,J), Z(IM+1,J) 
SRL Piva (151) 


AA: DROPLET HEIGHT; RR: DROPLET RADIUS. 


DIMENSION Z(IM+1,JM+1) ,R(IM+1,JM+1) 
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DIMENSION RSO(5*JS) 
COMMON /CONST/ ALFA, PI 
COMMON /TSTEP/ NT 

RRR=AA/RR 


DO 10 I=2,IM 
FI=1.-(I*1.-2.)/(IM*1.-2.) 
DO 10 J=2,JM 
GJ=(J*1.-2.)/(IM*1.-1.) 


IF(I.EQ.IM.OR.J.EQ.2) THEN 
IF (I.EQ. IM) THEN 

R(1,J) =RR*GJ 

Z(I,J)=0. 

ENDIF | 

IF(J.EQ.2) THEN 

R(I,J)=0. 

Z(1I,J) =AA*FI 

ENDIF 


ELSE 
FIR=1./ (FI*RRR) -FI*RRR 
A=(1.+GJ**2+2.*(1.-GJ) ) *RR/FIR 
B=-2.*(GJ+(1.-GJ)/GJ) /FIR 

IF (FIR.GT.1.E-3) THEN 

FACTOR=1. 

ELSE 

FACTOR=1.E-6 

ENDIF 

ALFA1=(1.+B**2) *FACTOR 

BEITA=2.* (A*B-(GJ+(1.-GJ) /GJ) *RR) *FACTOR 
GAMMA= (A**2+ (GJ**2+2.%*(1.-GJ) ) *RR**2) *FACTOR 
TF=BEITA**2-4. *ALFA1*GAMMA 

IF(TF.LT.0.)TF=0. 

IF (TF.GE.0.) THEN 
R(I,J)=(-BEITA-SQRT (TF) )/(2.*ALFA1) 

ELSE 


ENDIF 
Z(I,J)=SQRT((1./(FI*RRR)+FI*RRR) **2/4.*RR**2-R(I,J) **2) 
$ -(1./(FI*RRR) -FI*RRR) /2.*RR 

ENDIF 

10 CONTINUE 
DO 20 I=1,IM+1 
R(I,JM+1)=RR 
Z(1I,JM+1)=0. 

20 CONTINUE 

C..FOR IMAGINARY POINTS 

DO 30 I=1,IM+1 
R(I,1)=-R(I,3)+2-.*R(I,2) 
Z(I,1)=Z(I,3) 
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30 CONTINUE 
DO 40 J=1,JM 
R(IM+1,J7)=R(IM-1,7) 
Z(IM+1,J7)=-Z(IM-1,J)+2.*Z(IM,J) 

40 CONTINUE 
I=1 
AA1=AA+(Z(2,2)-Z(3,2)) 

=AA1/RR 
FI=1. 
Z(I,2)=FI*AAL 
DO 50 J=3,JM 
GJ=(J*1.-2.)/(IM*1.-1.) 
FIR=1./(FI*RRR) -FI*RRR 
A=(1.+GJ**2+2.%*(1.-GJ) ) *RR/FIR 
B=-2.*(GJ+(1.-GJ) /GJ) /FIR 
IF (FIR.GT.1.E-3) THEN 
FACTOR=1. 
ELSE 
FACTOR=1.E-6 
ENDIF 
ALFA1=(1.+B**2) *FACTOR 
BEITA=2.* (A*B-(GJ+(1.-GJ) /GJ) *RR) *FACTOR 
GAMMA= (A**2+ (GJ**2+2.*(1.-GJ) ) *RR**2) *FACTOR 
TF=BEITA**2-4.*ALFA1*GAMMA > 
IF(TF.LT.0.)TF=0. 
IF (TF.GE.0.) THEN 
R(I,J)=(-BEITA-SORT (TF) ) /(2.*ALFA1) 
ELSE 
PRINT *,'TF<0.??2??? IT IS WRONG! I=',I1,' J=',g 
ENDIF 
Z(I,J)=SQRT((1./(FI*RRR) +FI*RRR) **2/4.*RR**2-R(I,J) **2) 
$-(1./(FI*RRR) -FI*RRR) /2.*RR 


50 CONTINUE 
Carew s 
Rip -R (1,3) +2. *R( 1,2) 
Z(1,1)=2(1,3) 
C.5 2 Be ASSIGNING POINTS ON SOLID SURFACE WHERE.THE TEMPERATURES 
C....-ARE CALCULATED . 


IF (NT.EQ.1) THEN 
DO 60 J=1,JdS 
RSO(J)=.5*(R(IM,J+1)+R(IM,J+2) ) 

60 CONTINUE 

C...USE FINER GRID NEAR THE EDGE 

DO 61 J=JS+1,3*JS 
AJ=(1.*3-1.%*3S) /(2.*JS) *PI/2. 
RSO (J) =RR*(1.+2.*(1.-COS(AJ))) 

61 CONTINUE 
DR1=RS0 (3*JS) -RSO(3*JS-1) 
DO 63 J=3*JS+1,4*JS 
RSO(J)=RSO(J-1)+2.*DR1 

63 CONTINUE 
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DO 64 J=44IS+1,5%*IS 
RSO(J)=RSO(J-1)+8.*DR1 
64 CONTINUE 


ENDIF 
ais us sleie os 

RETURN 

END 
Cc 2 

SUBROUTINE DZ0(IM,JM,R,Z,RO3,Z03,DZ) 
Cc 


C....SUBROUTINE FOR FINDING THE DZ AT THE LIQUID AND SOLID 
C....INTERFACE 
DIMENSION R(IM+1,JM+1) ,Z(IM+1,JM+1) , DZ (JM) 
DIMENSION RO3 (IM) ,Z03 (IM) 
I=IM 
DO 135 JJ=2,IM 
IF (JJ.EQ.JM) THEN 
RM=RO3 (I) 
ZM=Z03 (I) 
=RO3 (I-1) 

ZP=Z03 (I-1) 
ELSE 
RM=.25*(R(I,JJ)+R(1,JJ+1)+R(I+1,JJ+1)+R(I+1,J5) ) 
ZM=.25*(Z(1,JJ)+Z(L,JI+1)+Z(I+1,J3J+1)+Z(I+1,JT)) 
RP=.25*(R(I-1,Jd)+R(I-1,JJ+1)+R(I,JJ+1)+R(1I,dJ) ) 
ZP=.25*(Z(I-1,3I)+Z (I-1,JI+1) +Z(I,JI+1)+Z(I,I)) 
ENDIF 
DZ (JJ) =SQRT( (RM-RP) **2+(ZM-ZP) **2) 

135 CONTINUE 
RETURN 
END 


SUBROUTINE SOLVE(N,A,B,X, IWK, WK) 


Clean hs SUBROUTINE FOR SOLVING MATRIX 
DIMENSION A(N,N),B(N) ,X(N) 
DIMENSION IWK(N) , WK(N+N*(N-1) /2) 
INTEGER N 
COMMON /WORKSP/RWKSP 
CALL L2NRG(N,A,N,A,N,WK, IWK) 

DO 100 J=1,N 

X(J)=0.0 

DO 100 I=1,N 

X(J)=X(J)+A(J,1I) *B(I) 
100 CONTINUE 


RETURN 
END 
Cc 
SUBROUTINE SURF(DTIME,IM,JM,RR,VO,AA,R,Z,DVDT, 
$ EWIS, COEFFH,SX,RHOW,CA) . 
G 
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een. SUBROUTINE FOR CALCULATING THE SIZE OF THE LIQUID DROP RR, AA 
con. WHILE RR IS A CONSTANT 
DIMENSION R(1+IM,1+JM) ,Z(1+IM,1+JM) , EWIS (JM) ,SX (JM) 
COMMON /MATER/ LMAT 
PI=4.*ATAN(1.) 
C....VO=PI*AA/6.*(3.*RR**2+AA*%2) 
SSX=0. 
DO 10 J=2,3M 
RJ=.5%*(R(2,J)+R(2,J+1) ) 
DRJ=SQRT ((R(2,J+1)-R(2,d) ) **2+(Z(2,d+1)-Z(2,70) ) **2) 
SSX=SSX+EWIS (J) *COEFFH*SX (J) *RJ*DRI 
10 CONTINUE 
DVDT=2.*PI*.624/ (RHOW*CA) *SSX 
VO=VO-DVDT*DTIME 
C...SOLVING CUBE EQ. AA**3+3.*RR*82*AA-6./PI*VO=0. 
U=EXP(1./3.*LOG(3.*VO/PI+SQRT (RR**6+9./PI**2*VO**2) )) 
- AA=U-RR**2/U 


SUBROUTINE NEXTR(NT,KT,FTIM,R, TIME,W,F,WF,RSO,W1) 
DIMENSION W(NP,NP) ,F(NP) , FTIM(NP, KT) ,W1(NP,NP, KT) 
DIMENSION WF(NP) ,R(NP) ,RSO(NP) 
COMMON /CONST/ ALFA, PI 
COMMON /GRID/ DRO,DTO,NDP,NP 
COMMON /INTG/ ATO,ERRREL, IRULE 
COMMON /MATER/ LMAT 
DO 149 II=1,NP 
R(II)=0. 
149 CONTINUE 

IF (NT.LE.KT) THEN 
KKT=NT 
ELSE 
KKT=KT 
ENDIF 
DO 150 IT=1,KKT 

C SELECT PROPER WEIGHTS 
TO=IT*DTO 

C..GET W FROM THE PREVIOUS CALCULATION 
IF (IT.LE.KKT-1.OR.KKT.EQ.KT) THEN 
DO 151 I=1,NP 
DO 151 J=1,NP 
W(I,J)=W1(1I,J,1T) 

151 CONTINUE 

ELSE 

C...CAL. THE WEIGHT FOR CURRENT TIME STEP 
CALL WEIGHT (TO,W,RSO) 

C...SAVE FOR LATER STEPS 
DO 152 I=1,NP 
DO 152 J=1,NP 
WIL (IT) HW; T) 
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152 CONTINUE 
ENDIF 
C PICK UP CORRESPONDING FLUX 
DO 100 I=1,NP 
F(I)=FTIM(I,IT) 
100 CONTINUE 
C SUMMATION IN SPACE 
CALL MULT(NP,NP,W,F, WF) 
C SUMMATION IN TIME 
CALL ADD(NP,WF,R) 
150 CONTINUE 
RETURN 
END 


SUBROUTINE MULT(M,N,W,F, WF) 
REAL W(M,N),F(N) ,WF(N) 


DO 50 I=1,M 

WF(I)=0.0 

DO 50 K=1,N 

WF (I) =WF (I)+W(I,K) *F(K) 
50 CONTINUE 

RETURN 

END 


SUBROUTINE ADD(NP,WF,R) 
REAL WF(NP) ,R(NP) 


DO 50 I=1,NP 

R(I)=R(1I)+WF(I) 
50 CONTINUE 

RETURN 

END 


SUBROUTINE WEIGHT (TO,W,RSO) 
REAL W(NP,NP) ,RSO(NP) 

REAL QDAG, BSIOE,ERF 
EXTERNAL Fl 


COMMON /CONST/ ALFA, PI 

COMMON /GRID/ DRO,DTO,NDP,NP 
COMMON /INTG/ ATO,ERRREL, IRULE 
COMMON /RRO/ R,RO 

COMMON /MATER/ LMAT 


DO 50 I=1,NP 
DO 50 J=1,NP 
IF(J.EQ.1) THEN 
DRO=2.*RSO (J) 
ELSE 
IF(J.EQ.NP) THEN 
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DRO=DRO 
ELSE 
DRO=.5*(RSO(J+1) -RSO(J-1) ) 
ENDIF 
ENDIF 
=RS0 (I) 
RO=RS0 (J) 
C SET UP-.(W) FOR TO < DTO 
ERRABS=0.0 
IF (TO.LT.DTO) THEN 
W(I,J)=QDAG(F1,ATO,.5,ERRABS, ERRREL, IRULE, RESULT, ERREST) 
W(1I,J)=RESULT 
END IF 
C SET UP (W) FOR TO > DTO 
IF(TO.GE.DT0O) THEN 
ARG1=DRO/SQRT (16*ALFA*TO) 
ARG3=(R-RO) **2/ (4*ALFA*TO) 
ARG4=0.5*R*RO/ (ALFA*TO) 
IF(I.NE.J) THEN 
W(I,J)=(1.0/SQRT(4.0*PI*ALFA) ) *RO*(TO**(-1.5))* 
& BSIOE(ARG4) *EXP(-ARG3) *DRO*DTO 
ELSE 
W(I,J)=(R/TO) *BSIOE(ARG4) *ERF (ARG1) *DTO 
END IF 


50 CONTINUE 


FUNCTION F1(TO) 
REAL BSIOE,ERF 

COMMON /CONST/ ALFA,PI 
COMMON /GRID/ DRO,DTO,NP 
COMMON /RRO/ R,RO 


ARG1=DRO/SQRT(16*ALFA*TO) 

ARG3=(R-RO) **2/ (4*ALFA*TO) 

ARG4=0.5*R*RO/ (ALFA*TO) 

E=0.0 

IF (ARG3.LE.100.0) E=EXP (-ARG3) 

IF (R.NE.RO) THEN 

F1=(1.0/SQRT(4.0*PI*ALFA) ) *RO*(TO*#*(-1.5))% 

& BSIOE (ARG4) *E*DRO 

ELSE 

F1=(R/TO) *BSIOE(ARG4) *ERF (ARG1) 

END IF 

RETURN 

END 
e 3 
SUBROUTINE CTEMPO(IM,JM,JE,TREF,TWATER, TSOLID,CT,TTO) 
C...SUBROUTINE CALCULATING THE CONTACT TEMPERATURE TO BE USED AT 


OTT 


C...THE SOLID-LIQUID INTERFACE 
DIMENSION TWATER(IM*JM+JE) , CT (JM) 
COMMON /SOLID/ SK,SC,SR 
COMMON /WATER/ COEFKW, CPW, RHOW 


SQRCKW=SQRT (RHOW*CPW*COEFKW) 
SQRCKA=SOQRT (SR*SC*SK) 
ape ; 
DO 10 J=1,JM 
CT (J) =(SQRCKW* (TWATER ( (IM-2) *JM+J) +TREF) 
+SQRCKA*TSOLID) / (SQRCKW+SQRCKA) 
TT=( (SQRCKW+SQRCKA) *100.-SQRCKW* (TWATER ( (IM-2) *JM+J) 
+TREF) ) /SQRCKA 
C..DIMENSIONLESS TEMP. 
TTO=(TSOLID-TREF) / (TT-TREF) 
ely 
IF(CT(1) .GE.100.) THEN 
PRINT *,'THE CONTACT TEMP. IS GREATER THAN 100C!'! 
PRINT *,' *#** TRY AGAIN *#+! 
PRINT *,'ENTER ANY NUMBER TO CONTINUE OR CTRL-C TO STOP! 
READ (5,*)NCT 
ENDIF 
10 CONTINUE 
PRINT *,'CONTACT TEMPERATURE=',CT(1) 
PRINT *,'DIMENSIONLESS TEMP: (TS-TAMB) /(TREF-TAMB) =',TTO 
RETURN 
END 
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This report describes the research performed during the period July 1990 - July 1991 under 

a joint research program between the Mechanical Engineering Department of the University of 
Maryland and the Building and Fire Research Laboratory of the National Institute of Standards 
and Technology. The research is conducted by Graduate Research Assistants of the ME 
Department under the joint supervision of Dr. di Marzo ({UMCP) and Dr. Evans (CFR - NIST). 

This joint research program was initiated in January 1985. The long term objective of the 
study of droplet-solid interaction is to obtain information applicable to the extinguishment 
of fire throug a droplet array (e.g. spray). The solids of concern include low thermal 
conductivity materials, typical of fire applications. Several important results were obtained 
in the first years’ of research. In particular, the modelling of the boundary condition at the 
liquid-vapor interface (at the droplet exposed surface) was validated with the data collected 
for water droplets evaporating on an aluminum block (diMarzo 1986a, 1986b, 1988). 
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